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Abstract 

We show that the Balitsky-JIMWLK equations proposed to describe non-linear 
evolution in QCD at high energy fail to include the effects of fluctuations in the gluon 
number, and thus to correctly describe both the low density regime and the approach 
towards saturation. On the other hand, these fluctuations are correctly encoded (in 
the limit where the number of colors is large) in Mueller's color dipole picture, which 
however neglects saturation. By combining the dipole picture at low density with 
the JIMWLK evolution at high density, we construct a generalization of the Balitsky 
hierarchy which includes the particle number fluctuations, and thus the pomeron 
loops. After an additional coarse-graining in impact parameter space, this hierarchy 
is shown to reduce to a Langevin equation in the universality class of the stochastic 
Fisher-Kolmogorov-Petrovsky-Piscounov (sFKPP) equation. This equation implies 
that the non-linear effects in the evolution become important already in the high 
momentum regime where the average density is small, which signals the breakdown 
of the BFKL approximation. 
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1 Introduction 

Over the last decade, much progress has been realized towards under- 
standing the dynamics of QCD at high energies, in the vicinity of the unitar- 
ity limit. New theoretical approaches have been developed which encompass 
and extend the linear BFKL equation [1] originally proposed to describe high 
energy evolution in the leading-logarithmic approximation (with respect to 
Ins), as well as the earlier attempts to improve over the BFKL equation by 
including unitarity corrections (or gluon saturation) in the form of non-linear 
terms [2,3,4]. 

In the mid nineties, Al Mueller [5,6] has given an elegant construction 
of the BFKL wavefunction of an energetic hadron, the 'color dipole picture', 
which exploits the large-iV c approximation (with N c the number of colors) to 
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replace gluons by color dipoles as the effective degrees of freedom at small 
x, and gluon radiation by dipole splitting as the mechanism for evolution. 
The average number of dipoles obeys BFKL equation, but the overall descrip- 
tion goes beyond the original BFKL approach by including the many-body 
correlations generated through dipole splitting. Accordingly, the evolution is 
expressed as a hierarchy of equations for the dipole density correlations. The 
dipole picture cannot describe interactions among the dipoles in the wavefunc- 
tion (since such interactions would lead to color configurations with higher 
multipolar moments), so it cannot accommodate gluon saturation. But it can 
be used to study the onset of unitarity corrections in the elastic scattering be- 
tween two sets of dipoles which have evolved separately, provided the collision 
is viewed in the center-of-mass frame and the total energy is not too large 
[6]. (The unitarization is brought in by the multiple scattering of several pairs 
of dipoles from the two colliding systems.) More generally, the dipole picture 
can be used to describe the dilute (in the sense of non-saturated) part of a 
hadron wavefunction up to arbitrarily high energies, and this is how we shall 
actually use it throughout this paper. 

A different formalism to study unitarization at high energy, in which the 
unitarity corrections are encoded directly in the evolution equations, has been 
developed by Balitsky [7]. This involves a hierarchy of equations describing 
the evolution of the scattering amplitudes for the collision between a generic 
target and a set of relatively simple 'projectiles', which transform into each 
other under the evolution. The target is not evolving, but is generally taken to 
be a high-density system, or 'color glass condensate', represented by a strong 
color field. The evolution is achieved by boosting the projectile, but in such 
a way that the latter remains non-saturated, which implicitly restricts the 
maximal energy allowed. Thus, the evolution of the projectile wavefunction 
is in fact linear — for a dipole projectile and in the large-Ac limit, this is 
precisely the dipole picture discussed above 2 — , but Balitsky equations are 
nevertheless non-linear, since written for the scattering amplitudes : the non- 
linear effects correspond to multiple scattering between the components of the 
projectile and the color field in the target. 

Alternatively, one can keep the projectile unevolved, but use the increase in 
the energy in order to boost the target, and then study the small-x evolution 
of the strong color fields there. In doing so, some of the non-linear effects 
which were interpreted as multiple scattering from the perspective of projectile 
evolution will now appear as saturation effects in the evolution of the target. 
But to study this, one needs a formalism capable to deal with non-linear 
effects in the evolution of the wavefunction. For the reasons explained before, 



2 Recently, it has been explicitly shown, by Levin and Lublinsky [8], that the large- 
N c version of the Balitsky equations follows directly from the dipole picture for the 
projectile, after allowing for multiple interactions with the target. 
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this cannot be a large-iV c formalism in terms of dipoles, but rather it must 
allow for all the possible color multipoles, and the simplest and most natural 
way to do so is to work in terms of gluons. So far, the only formalism of this 
type is the color glass condensate (CGC) [9,10,11,12,13,14] (for reviews see 
Ref. [15]), in which the small-x gluons are described as classical color fields 
generated by color sources at larger values of x, which are randomly distributed 
(in color and coordinate space) with a weight function which evolves with the 
energy. This evolution is governed by a functional Fokker-Planck equation for 
the weight function, the JIMWLK equation [12,13,16], which is equivalent to 
an infinite hierarchy of equations for the correlation functions of the classical 
color fields. When applied to the scattering between a simple projectile and the 
CGC, the JIMWLK evolution reproduces Balitsky equations for the scattering 
amplitudes. In what follows, we shall refer to these equations as the Balitsky- 
JIMWLK equations. 

The Balitsky- JIMWLK equations will be further discussed in Sect. 3 be- 
low, but here we would like to emphasize a few important points: i) The 
evolution described by these equations is stochastic : Through the non-linear 
terms present in these equations, n-point functions with different values of n 
will mix under the evolution, thus generating new correlations with increasing 
energy, ii) The non-linear terms are associated with the presence of strong 
classical fields (the field of the target in the approach by Balitsky, or the field 
created by color sources at higher values of x in the CGC). Therefore, in 
the dilute regime where the fields are weak, the equations can be linearized, 
and once we do so, they decouple from each other (so that the various n-point 
functions obey independently the BFKL evolution). This should be contrasted 
to the dipole picture, where correlations beyond the BFKL equation appear 
already in the dilute regime. Hi) Even in the full equations with the non-linear 
terms included, the stochastic aspects turn out to become inessential at large 3 
iV c . For instance, in the absence of correlations in the initial conditions, the 
large-iV c version of the Balitsky hierarchy boils down to a single, non-linear, 
equation, which is deterministic. This is the equation originally derived by 
Kovchegov [17] and generally referred to as the Balitksy-Kovchegov (BK) 
equation. This suggests that the stochastic aspects of the Balitsky-JIMWLK 
equations are to be attributed to color fluctuations in the high-density regime, 
which are suppressed at large N c , unlike the colorless fluctuations in the par- 
ticle number, as encoded in the dipole picture, which are present already in 
the dilute regime and for large N c . 



3 By "large--/V c " we understand here the high-energy version of the large--/V c ap- 
proximation, due to Mueller [5,6], in which gluon exchanges which are suppressed 
by factors of are still counted as leading order effects provided they are en- 

hanced by appropriate powers of the energy. The neglected terms are those which, 
for a given power of the energy, are suppressed by higher powers of 1/N% than the 
leading-order terms. 
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The previous discussion — which will be further substantiated and illus- 
trated with a few Feynman diagrams in Sect. 3 — points out towards an 
insufficiency of the Balitsky-JIMWLK equations, which fail to include the 
correlations associated with fluctuations in the particle number. Because of 
that, the mean field approximation (the BK equation) should work very well 
for these equations at large iV c and for uncorrelated initial conditions (like 
the scattering off a large nucleus), and this is indeed what is seen in numeri- 
cal studies of JIMWLK evolution [18]. So far, most studies of saturation and 
unitarity have in fact focused on the BK equation, which is much easier to 
handle in practice than the general equations, and whose properties are by 
now rather well understood [19,20,21,22,23,24,25,26,27,28,29]. 

However, there were indications already in the dipole picture (see espe- 
cially the numerical studies by Salam [30,31]) that the effects of fluctuations 
should be important, in particular, in the evolution from a low density regime 
to a high density one, and thus in the approach towards saturation. Indeed, 
in the dilute regime, the particle number fluctuations are a natural mecha- 
nism to generate higher-point density correlations, which then evolve into the 
non-linear terms responsible for saturation. E.g., if one starts with an isolated 
gluon (or dipole) at high transverse momentum (k±), as radiated from the bulk 
of particles at lower k±, then higher-point correlations at high k± — which 
were originally absent — will get built in the early stages of the evolution, 
as correlations in the splitting of the original dipole or its close descendence. 
Later on, these correlations get amplified by the standard BFKL evolution, 
and eventually influence the evolution of the lower-point correlations when 
the density is high enough. 

The importance of fluctuations in the evolution towards saturation has 
been reiterated in recent publications [32,33,34]. In Ref. [32] it has been ar- 
gued that the fluctuations slow down the approach towards the unitarity limit 
as compared to the mean field approximation (MFA) [14,19]. In Ref. [33], a 
modified MFA has been proposed, in which the unitarity constraints associated 
with fluctuations have been simulated by imposing a barrier at high momenta 
on the BFKL equation. The solution to the ensuing equation shows that fluc- 
tuations reduce considerably the rate of growth of the saturation momentum 
with the energy (or 'saturation exponent'), and that the corresponding mean- 
field limit is reached only very slowly, logarithmically, as a s -> 0. Finally, in 
Ref. [34], the event-by-event picture of the evolution has been considered, 
in which fluctuations appear naturally because, in a given event, the parti- 
cle occupation numbers are discrete. Based on this picture, a correspondence 
has been established between high-energy evolution in QCD and a class of 
stochastic particle models which are actively studied in statistical physics. By 
using known results for the latter, the authors of Ref. [34] have confirmed the 
result in Ref. [33] about the slow convergence of the saturation exponent to its 
mean-field value, and further argued that fluctuations should eventually wash 
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out the geometric scaling [35] property of the average amplitude, which was 
known to hold for the BK equation [26,27,28,29]. However, the arguments in 
Ref. [34] cannot predict how fast (i.e., after what energy evolution) is geometric 
scaling violated if one starts with an initial condition that shows scaling, nor 
they allow one to study the preasymptotic evolution at intermediate energies 
and for realistic values of a s . 

More importantly, there was some confusion in Ref. [34] about what should 
be the right formalism to study this evolution. Namely, it has been assumed 
there that the correct evolution law is the one provided by the Balitsky- 
JIMWLK equations, but this cannot be right, since these equations miss the 
relevant fluctuations, as we have argued above. It is then natural to ask, what 
is/are the equation(s) that one should solve (e.g., numerically) in order to 
study the interplay between fluctuations and saturation, and thus be able to 
verify the results in Refs. [33,34] and improve over them. In other terms, what 
is the correct formalism to describe evolution in QCD at high energy ? 

It is the main purpose of this paper to provide an answer to this question, 
in the form of a new set of evolution equations which generalize the (large- 
iV c version of the) Balitsky-JIMWLK equations by approximately including 
the effects of particle number fluctuations. Our basic observation is that, in 
fact, we do have the necessary ingredients to describe both fluctuations and 
saturation — these are the color dipole picture, and the JIMWLK evolution, 
respectively — and that, although these ingredients may look too different 
from each other to be simply married in a unified theoretical description of 
the lightcone wavefunction, they can still be naturally merged with each other 
in the evolution equations for scattering amplitudes. Since relying on the dipole 
picture, our subsequent construction is limited to the large-iV c limit, to which 
we shall restrict ourselves throughout the following analysis. It remains as 
an open problem at this stage whether it would be possible to develop a 
wavefunction formalism which incorporates all these physical ingredients (e.g., 
through a suitable extension of the JIMWLK formalism), and thus rederive 
the results for scattering amplitudes that we shall present here, together with 
their generalization to arbitrary N c . 

Let us briefly anticipate here the main steps in these construction, and 
explain the other approximations that we shall need. After briefly explaining 
the physical picture of particle number fluctuations in the dipole picture in 
Sect. 2, we shall then construct, in Sect. 5, the equations which describe the 
evolution of the many-body density correlations in this formalism. To that 
aim, we shall use the recent formulation of the dipole picture in Refs. [36,37], 
which is more convenient for a study of fluctuations since it follows explicitly 
the evolution of a given configuration of N dipoles. The equations that we 
shall obtain in Sect. 5 have been already presented by Levin and Lublinsky in 
a very recent publication [8], but our respective derivation will be somewhat 
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different, and also better suited for our present purposes. 

In the meantime, in Sect. 4, we shall also analyze a simple stochastic 
particle model borrowed from statistical mechanics (a zero-dimensional ver- 
sion of the model in Ref. [38]), which in spite of its formal simplicity has the 
conceptual advantage over the dipole picture to include both fluctuations and 
recombination (and thus saturation). The explicit manipulations permitted by 
the simplicity of this model will help us developping some physical intuition, 
and will also serve as a guidance for the corresponding manipulations in QCD, 
where such a unified theoretical description of the wavefunction is still lacking. 
A more elaborated, one-dimensional, version of the same model (the one which 
is actually considered in Ref. [38]) will be then discussed in the Appendix. This 
will give us the opportunity to introduce the stochastic Fisher-Kolmogorov- 
Petrovsky-Piscounov (sFKPP) equation, which in this context emerges as an 
exact equation (mathematically equivalent to the original model) , and which 
will also emerge as an approximation to the corresponding equations in QCD. 

In Sect. 6 we construct the new equations for scattering amplitudes in 
QCD, which are the main results of this paper. First, in Sect. 6.1, we relate 
dipole correlations in the target to scattering amplitudes for projectile dipoles, 
by using approximations which make sense in the dilute regime. This allows 
us to translate, in Sect. 6.2, the equations for densities constructed in Sect. 5 
into corresponding equations for the scattering amplitudes. These equations 
are still linear, as appropriate in the low-density regime, but they include 
the effects of dipole number fluctuations in the target (so they form already a 
non-trivial hierarchy). Then, these equations are extended to the high-density 
regime by adding the non-linear terms expected from the JIMWLK evolution 
(i.e., the same non-linear terms as in the Balitsky hierarchy). Because of the 
combined effect of fluctuations ('gluon splitting') and recombination ('gluon 
merging'), the final equations generate pomeron loops through iterations. In 
Sect. 6.3, we show that, after a coarse-graining in the impact parameter space, 
the whole hierarchy can be equivalently replaced by a specific Langevin equa- 
tion, which is formally the BK equation supplemented with a multiplicative 
noise term. 

In our last section 7, we explore some physical and mathematical conse- 
quences of the new equations. In Sect. 7.1, we show that, mediating a gra- 
dient expansion of the BFKL kernel known as the 'diffusion approximation', 
the Langevin equation of Sect. 6 can be cast into the form of the sFKPP 
equation. The latter emerges as an effective equation for a variety of prob- 
lems in physics, biology, and chemistry (see Refs. [39,40] for recent reviews), 
and it has been extensively studied in recent years in the statistical physics 
literature. By using some known results about this equation, we confirm (in 
Sect. 7.2) previous findings in Refs. [33,34], and thus conclude that the present 
Langevin equation is indeed the evolution law underlying the physical discus- 
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sion in Ref. [34]. Finally, in Sect. 7.3, we discuss a rather dramatic consequence 
of the stochastic nature of the evolution, which entails the breakdown of the 
BFKL approximation in the high momentum regime, where the average gluon 
density, or the average scattering amplitude, are small. Because of fluctua- 
tions, the saturation momentum becomes a random quantity which can take 
different values from one event to another, or from one impact parameter 
to another (for the same event). Then the evolution develops high-density 
'spots' where the gluons are at saturation even for relatively high transverse 
momenta, well above the average saturation momentum. At sufficiently high 
energy, all the correlations are dominated by such dense spots, and the mean 
field approximation (i.e., BFKL equation at high k±) breaks down. This may 
look unexpected, but in fact it can be related to a similar behavior observed 
by Salam in his Monte-Carlo studies of the dipole picture [30,31]. The 'dense 
spots' were clearly seen in those numerical simulations, but of course in that 
context the local density was never saturating. 



2 Physical motivation 

Given a high energy hadron with rapidity Y, we would like to under- 
stand fluctuations in the tail of the gluon distribution at transverse momenta 
k± ^> Q S (Y) (with Q S (Y) the saturation momentum), and, in particular, the 
influence that such fluctuations may have — through their subsequent evolu- 
tion — on the approach towards saturation (or towards unitarity in the colli- 
sion with an external projectile). The high-Zc^ tail is a priori a dilute regime in 
which the gluon occupation numbers are small, of 0(1), so we expect impor- 
tant fluctuations associated with the discreteness of the particle (here, gluon) 
number. Given that the gluon number density is represented, in a field the- 
oretical formulation, by the two-point correlation function of the color fields 
in the hadron wavefunction, we deduce that the relevant fluctuations should 
be encoded in the four-point, and higher, correlations. Alternatively, and sim- 
pler, in the large-iV c limit to which we shall stick throughout this paper, the 
tail of the distribution can be described in Mueller's dipole picture [5,6], in 
which gluons are effectively replaced by quark-antiquark pairs of zero trans- 
verse size and in a color octet state. The 'color dipoles' then emerge as color 
singlet states built with the quark component of some gluon and the antiquark 
component of some other gluon with roughly the same rapidity. In this pic- 
ture, which applies so long as the gluon density is low enough for saturation 
effects to be negligible, the (unintegrated) gluon distribution is measured by 
the average dipole number density (n(r,Y)) which obeys BFKL equation [1]. 
But the dipole picture goes beyond the strict BFKL equation by including cor- 
relations in the distribution of dipoles, namely those correlations which follow 
from dipole splitting in the course of the evolution. The simplest such cor- 
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relations, which also encompass the particle number fluctuations that we are 
mainly interested in, are those encoded in the pairwise dipole number density 
(n( 2 )(n, r2, y)) which is a measure of the probability to find simultaneously 
two dipoles, with transverse sizes r\ and r 2 , respectively (see Sect. 5 below for 
a more precise definition). 

Since it is preferable to consider quantities which are measurable (at least, 
in principle), we shall study the effect of fluctuations on the scattering be- 
tween the dilute hadronic target (described within the dipole picture) and an 
external projectile, which is itself chosen as a set of dipoles in some fixed con- 
figuration. This choice is convenient since, as well known, the dipole-dipole 
scattering is quasi-local in phase-space: a dipole projectile essentially counts 
the numbers of dipoles in the target having the same transverse size 4 and 
impact parameter as itself. Specifically, the scattering amplitude for a single 
dipole can be estimated as T(r,b,Y) ~ a 2 f(r,b,Y), where r and b denote 
respectively the size and the impact parameter of the projectile dipole, a 2 
measures the scattering amplitude for two dipoles with roughly the same size 
and nearby impact parameters, and f(r,b,Y) is the (dimensionless) dipole 
occupation factor in the target, and is related to the corresponding number 
density via /(r, b, Y) ~ r 4 n(r, b, Y) (see Sect. 6). 

To be sensitive to fluctuations, the projectile must involve at least two 
dipoles, and here we shall consider the case where it contains exactly two. 
Both dipoles are assumed to be small, so the individual scattering amplitudes 
(Ti) ~ a 2 s (fi) (with (Ti) = (T(r h b h Y)) and i = 1, 2) are both small: (Tj) < 1. 
The question is then, how large can be (T^) ~ a 4 (/i/ 2 ) (the scattering 
amplitude for the simultaneous scattering of both dipoles) ? Clearly, in the 
absence of correlations in the target, or if the correlations are only weak, one 
has (/i/ 2 ) ~ (/i)(/ 2 ), and then (T^) ~ (T\)(T 2 ) is much smaller than either 
(Ti) or (T 2 ). However, as we discuss now, there are kinematical situations 
in which one expects strong correlations among the dipoles in the target, 
which will drastically enhance (T^) with respect to its estimate (Ti)(T 2 ) for 
incoherent scattering. 

The first such a situation is not really relevant for the main stream of this 
paper, and the only reason for mentioning it here is to distinguish it from the 
more interesting case to be discussed after. Namely, this is the situation in 
which both external dipoles scatter off the same dipole in the target. This is 
advantageous in the very dilute regime where (/) <C 1 (or (T) <C a 2 ), since for 
such a process there is a single low-density penalty factor, and the correspond- 
ing contribution (T^) ~ aj(f) is indeed much larger than the uncorrelated 



4 More precisely, (quasi)locality in the dipole size holds so long as the dipole distri- 
bution in the target is in a genuine BFKL regime, i.e., it is characterized by some 
'anomalous dimension' ; see Sect. 6 for details. 
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piece ~ aj{f) 2 . On the other hand, in view of the previous considerations 
on the quasi-locality of dipole-dipole scattering, it is clear that this situation 
requires a rather fine tuning between both the transverse sizes and the im- 
pact parameters of the two incoming dipoles (r 1 ~ r 2 and b\ ~ b 2 ), which as 
we shall see is not an interesting configuration for the evolution of the gluon 
distribution to high energies. 

The second situation, which is more relevant for our purposes here, is the 
one in which the two dipoles making the projectile are contiguous in transverse 
space, meaning that the quark leg of one dipole is close to the antiquark leg of 
the other. If we use X\ and y\ [x 2 and y 2 ) to denote the transverse positions 
of the quark and respectively antiquark leg of the first (second) dipole, then 
contiguous configurations correspond to either yi as x 2 or to X\ ~ y 2 . Such 
configurations are interesting because the scattering amplitude (T^) for two 
contiguous dipoles enters the evolution equation for the scattering amplitude 
of a single dipole in the regime where unitarity corrections become important 
(see Sect. 3). The reason why (T^) is enhanced as compared to (T) 2 for 
such configurations has to do with the dynamics of the evolution in the dipole 
picture: When increasing rapidity from Y to Y+dY , a dipole with legs at x and 
y can evolve by radiating one soft (i.e., small-x) gluon located at z, a process 
which at large- N c is tantamount to the original dipole (x, y) splitting into two 
new dipoles (x, z) and (z, y), which are contiguous. This evolution leads to 
an increase in the scattering amplitude of a projectile made of two contiguous 
dipoles (with appropriate transverse coordinates) which is proportional to the 
average density (n(x, y)) Y of the parent dipole at rapidity Y. That is, it is 
the dipole density (n), rather than the dipole pair density (n^), which acts 
directly as a source for (T^), and this is certainly advantageous in the dilute 
regime where 5 (n^) <C (n). 

For instance, if one starts with a single dipole (x , y ) at Y — 0, then the 
only way to find a pair of dipoles after an evolution dY is that this pair be made 
of the dipoles (x , z) and (z, y ), with arbitrary z. More generally, the Y- 
evolution of an arbitrary target will generate "high-/cj_ fluctuations" (i.e., small 
dipoles which split off the larger, preexisting dipoles) with low occupation 
numbers, of 0(1), and the only way to increase the dipole occupancy in those 
high-Zc^ bins in the next few steps of the evolution is through the splitting 
of the original fluctuations. It is only when the average pair density (n^) 
becomes large enough that the normal BFKL evolution (here, for (n^)) takes 
over, and the fluctuations (at that particular value of k±) cease to play a role. 
We see that a high-/cj_ fluctuation play the crucial role of a seed for extending 
the gluon distribution towards larger transverse momenta with increasing Y. 



5 More correctly, this inequality should be written for the dimensionless occupation 
numbers, as (//) < (/). 
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Fig. 1. The geometry of dipole splitting 



It is in fact easy to turn some of the above considerations into more explicit 
formulae: The differential probability for a dipole (x, y) to split into new 
dipoles (x, z) and (z, y) is known within perturbative QCD (to leading-log 
accuracy with respect to Y — Inl/x) as (a s = a s N c /ir) 



a s (x-y) 2 



dP(x,y\z) = -± V y} <PzdY. (2.1) 
2n (x — z) 2 (z — yy 

This formula, together with elementary geometry considerations (see Fig. 1), 
immediately imply (with v\ = X\ — y\ and r 2 = X2 — y 2 ) : 



d i (2), \\ a s (r-i + r 2 ^ 2 



{S i2) (x 2 - yi)(n(aji, y 2 ))y 



fiuct. 27r r\r\ 

+ 5 ( - 2 \y 2 -x 1 )(n(x 2 , Vl )) Y ) .{2.2) 

As indicated in the l.h.s. of Eq. (2.2), this is just the "fluctuating" contribution 
to the evolution of (n^), in which the two measured dipoles are generated via 
the splitting of a unique original dipole. The general equation for d(n^) / 8Y , 
to be constructed in Sect. 5 within the dipole picture, involves also terms 
linear in (n^), which correspond to the usual BFKL evolution. By also us- 
ing the relation T(r, b, Y) ~ a 2 r 4 n(r, b, Y), between (target) dipole densities 
and (projectile) dipole scattering amplitudes, it is possible to transform evo- 
lution equations for the dipole densities, so like Eq. (2.2), into corresponding 
equations for the scattering amplitudes. This will be discussed in Sect. 6. 



3 The Balitsky— JIMWLK equations and beyond 



Our discussion of the Balitsky-JIMWLK equations in this section will be 
rather schematic, as our intention is merely to show — by inspection of the 
structure of these equations, and with the help of a few Feynman diagrams - 



that these equations do not include the effects of gluon number fluctuations 
in the target wavefunction. Since, as we shall later argue, these fluctuations 
serve as the initiators of the "pomeron loops" in the target, our conclusion 
also implies that the evolution generated by these equations fails to include 
the pomeron loops. This failure was already recognized in the literature in 
relation with the Kovchegov equation, which is a mean field approximation 
to the Balitsky-JIMWLK equations, but to our knowledge it has not been 
argued so far for the general equations. 

For consistency for the other developments in this paper, we shall restrict 
also the present discussion to the (high-energy version of the) large- N c limit, 
in which Balitsky-JIMWLK equations close in the space of dipoles. That is, 
they reduce to a hierarchy of equations for the evolution of the scattering am- 
plitudes of a set of N dipoles, with N — 1, 2, . . . , which scatter off a generic 
target. The color fields in the target can be strong — that is, the target can 
be at saturation (a 'color glass condensate') — , but the equations include 
non-linear effects which ensure that scattering is unitary. As we shall see, the 
non-linear effects in the evolution can be interpreted as either multiple scat- 
tering, or saturation, depending upon the perspective from which one views 
the evolution (as projectile or, respectively, target evolution). 

The simplest way to present the (dipolar version of the) Balitsky-JIMWLK 
equations is to notice that the whole hierarchy can be generated from the fol- 
lowing "operator" equation: 



dT Y (x,y) _a s r 2 (x-y) 2 

dY ~2^J dZ {x-zf{y-zf [6A) 

{ - T Y (x, y) + T Y (x, z) + T Y (z, y) - T Y (x, z)T Y (z, y)}. 

By "operator equation" we simply mean that, in order to deduce the equa- 
tion satisfied by the iV-point function (T (7V) ) y = (T(1)T(2) • • -T(N)) Y (with 
T{i) = T(xi,yi)), it is sufficient to multiply Eq. (3.1) by appropriate powers 
of T and then use Leibniz' rule; e.g. : 



£<™)> y = <^r(2)) y + <r(l)^) y . (3.2) 

In particular, the equation obeyed by the average scattering amplitude of a 
single dipole is immediately obtained as: 



>^5/ A |,y-,y (3 - 3) 

( - T(x, y) + T(x, z) + T(z, y) - T(x, z)T(z, y)) y . 



12 



As anticipated in Sect. 2, the r.h.s. of this equation involves the scattering 
amplitude (T^) = (T(x, z)T(z,y)) for two conti guous dipoles. More gener- 
ally, the equation for (T^) Y involves also (T^ N+1 ^) Y , so Eq. (3.1) generates 
an infinite hierarchy of equations which decouple from each other only in the 
weak scattering regime, where (T^ N+1 ^) <C (T^') and each iV-point function 
obeys separately the BFKL equation for evolution in any of its N arguments. 

But the hierarchy generated by Eq. (3.1) is sufficiently simple for the 
corresponding evolution to be quasi-deterministic. It is indeed easy to check 
that, if the initial conditions at Y — Y are chosen in factorized form, i.e., 
(T(l) • • -T(N))q = (T(l))o ■ • • (T(N))o, then this factorized form is preserved 
by the evolution up to arbitrarily large Y : schematically, (T^)y = (T) Y 
with the one-point function (T) Y obeying the BK equation (i.e., the equa- 
tion obtained by replacing T — > (T) Y into Eq. (3.1)). More generally, it has 
been shown in Refs. [41,42] that the hierarchy generated by Eq. (3.1) admits 
a one-parameter family of fully factorized exact solutions. This strongly sug- 
gests that, in their simplified form valid at large N c , the Balitsky-JIMWLK 
equations do not generate new correlations, but only propagate those already 
encoded in the initial conditions. This simplifying feature does not hold also 
for the full equations, which include additional, multipolar, operators. But the 
correlations induced by these operators are suppressed by powers of 1/N%, and 
thus cannot be associated with fluctuations in the gluon (or dipole) number. 
Rather, as discussed in the Introduction, they describe color fluctuations. 

This conclusion is further substantiated by an analysis of the diagrammatic 
content of the Balitsky-JIMWLK equations, to which we now turn. At this 
point one should recall that the structure of the perturbation theory, and thus 
the form of the diagrams, depends upon the frame in which one is viewing the 
evolution: 

a) In the original derivation by Balitsky [7], the evolution is implemented 
by boosting the projectile, which then evolves through (small-rr) gluon radia- 
tion. In the large-iV c limit in which the projectile is a collection of N dipoles, 
its evolution amounts to the splitting of one of these dipoles into two new 
dipoles, followed by the interaction between the final system of N + 1 dipoles 
and the target. Then Eq. (3.1) applies to one such a dipole which has split, 
and the various terms there are easily interpreted: The two linear terms with 
positive sign, T Y (x,z) and T Y (z,y), describe the independent scattering of 
the daughter dipoles with the target, the quadratic term with a negative sign 
corrects for an overcounting of their simultaneous scattering, and the linear 
term with a negative sign is the "virtual term" which expresses the possibility 
that the parent dipole (x, y) survive without splitting. 

In this picture, the non-linear terms in the evolution are thus associated 
with multiple scattering, but they indirectly reflect saturation effects in the 
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(a) (b) 




(e) (f) 

Fig. 2. Diagrams for single dipole scattering: (a) the tree-level contribution; (b,c,d) 
one step evolution of the projectile; (e,f) one step evolution of the target. 

target. But as far as the projectile is concerned, there is still no saturation: 
The N dipoles composing the projectile are not allowed to interact with each 
other, which is a good approximation only so long as a 2 s N <C 1. Since the 
total number of dipoles within the projectile grows exponentially with Y, this 
is clearly a low energy approximation. 

In Fig. 2, we display a few Feynman diagrams which illustrate the processes 
encompassed by Eq. (3.1) from the perspective of projectile evolution. These 
diagrams are proper to the dipole picture, so the evolution is shown as dipole 
splitting rather than as gluon emission. For more clarity, we show only one 
diagram contributing to each type of process, which is moreover taken at the 
lowest non-trivial order in perturbation theory. Thus, the scattering between 
a dipole and the CGC target starts at two gluon exchange, as shown in Fig. 



14 



small-x gluon 



Fig. 3. A typical gluon cascade which contributes to a classical color held conhgu- 
ration in the color glass condensate 

2. a. The one-step evolution of this amplitude then involves the diagrams in 
Figs. 2.b, c, d : Fig. 2.b is the virtual term (the original dipole interacts 
before splitting), Fig. 2.c shows the scattering of one of the daughter dipoles 
(there is a similar diagram for the other dipole), while Fig. 2.d describes their 
simultaneous scattering. We shall return to a discussion of these diagrams 
after also introducing the picture of the target evolution. 

b) A different derivation of the Balitsky equations is provided [13] by the 
color glass condensate (CGC) formalism, in which the evolution is achieved 
by boosting the target, and the non-linear terms correspond explicitly to 
saturation effects in the target wavefunction. The 'color glass' is a random 
superposition of classical color fields with a gauge-invariant weight function 
which evolves with Y according to JIMWLK equation [12,13,16]. Diagramat- 
ically, each classical field configuration corresponds to a set of gluon cascades 
which are initiated by 'color sources' (e.g., valence quarks) carrying a sizeable 
fraction of the target rapidity Y, and which end up with a small-x gluon at 
the comparatively low rapidity of the projectile. (See Fig. 3 for a pictorial 
representation.) When increasing Y, these cascades evolve through gluon ra- 
diation from the classical field created in the previous steps. For weak fields 
(A & 1), corresponding to low gluon density, the JIMWLK equation reduces 
to the standard BFKL evolution: the various iV-point functions evolve inde- 
pendently from each other, and rise rapidly with Y. But at high density, the 
classical fields are strong (A ^> 1, and even A ~ 1/g in the 'condensate' regime 
at saturation), and the dynamics is fully non-linear. 

For the subsequent discussion, it is important to notice that the non- 
linearities in the CGC correspond to gluon merging, but not also to gluon 
splitting : As illustrated in Fig. 3, different gluon cascades can fuse with each 
other, via the non-linear effects encoded in the classical field equations (the 
Yang-Mills equations), and also in the emission and the propagation of the 
quantum gluons. This recombination process is what we refer to as 'gluon 
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merging'. On the other hand, the quantum gluons radiated from the classical 
fields are not allowed to further radiate small-x gluons by themselves; that is, 
the 'gluon splitting' is not included in the evolution. This amounts to neglect- 
ing gluon number fluctuations in the CGC picture. (Indeed, gluon splitting 
would be the CGC analog of the dipole splitting discussed in Sect. 2.) This 
is a good approximation at high density, where gluon splitting is indeed sup- 
pressed as compared to direct emission from the strong classical fields. But as 
we shall argue later, gluon splitting is a dominant effect in the dilute regime at 
high k±, and plays also an essential role in the approach towards saturation. 

This discussion has an obvious consequence for the structure of non- 
linearities in the JIMWLK evolution equations for Green's functions: gluon 
merging can only reduce the number of active fields in a correlation function. 
Therefore, the r.h.s. of the equation for the iV-point function involves in gen- 
eral all the n-point functions with n > N, but not also those with n < N 
(compare in this respect with Eq. (2.2), where the opposite situation happens). 

Consider now the scattering between the CGC and an external projectile 
which is a collection of dipoles. The scattering amplitude for a single dipole 
can be computed in the eikonal approximation as 



T(x,y) = 1 - w tr(v^x)V(y)), (3.4) 

where V\x) and V(y) are Wilson lines describing the scattering of the quark 
and, respectively, antiquark component of the dipole, e.g. (with P denoting 
path-ordering in a; - ) : 

V ] {x) = Pexp (ig J dx~A+(x-,x)t a ^ , (3.5) 

and is the longitudinal component of the color field in the target (the only 
dynamical field left in the problem), whose correlations are described by the 
CGC. The evolution equation for (T(1)T(2) • • -T(N)) Y is then obtained by 
using the definition (3.4)-(3.5) of T together with the equations satisfied by the 
field correlators (A + (l)A + (2) ■ ■ -)y, which in turn follow from the JIMWLK 
equation 6 . As anticipated, the ensuing equations are precisely those originally 
obtained by Balitsky, and which at large N c are generated by Eq. (3.1). 

Let us now describe the diagrammatic interpretation of these equations 
from the perspective of target evolution. Once again, we show only a minimal 
set of representative diagrams. To easy read these diagrams, notice that the 

6 In practice, it is often simpler to work out directly the evolution equations for 
the correlations of the Wilson lines. But thinking in terms of field correlations helps 
clarifying the diagrammatic interpretation of the JIMWLK evolution. 
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lowest-order contribution to the scattering amplitude (3.4), as obtained after 
expanding the Wilson lines there, is the 2-gluon exchange: 

2 

T(x,y) « ^(A+(x)-A+(v)) 2 + 0(g 3 ). (3.6) 

The two diagrams shown in Figs. 2.e and f correspond to the linear and non- 
linear terms in Eq. (3.1), respectively. The diagram in Fig. 2.e represents the 
first step in the BFKL evolution of the 2-point function (A + A + )y (the 'unin- 
tegrated gluon distribution'); in terms of projectile evolution, it corresponds to 
the diagrams shown in Figs. 2.b and c. But Fig. 2.f exhibits a non-linear effect 
which goes beyond BFKL evolution: four gluons merge into two (a 4-point 
function reduces to a 2-point function) through a vertex which at large-iV c 
can be recognized 7 as the 'triple pomeron vertex'. Fig. 2.f corresponds to Fig. 
2.d for projectile evolution; both diagrams involve the same triple pomeron 
vertex, only with a different interpretation: i) In Fig. 2.d, this vertex describes 
the splitting of one dipole into two dipoles which then scatter both with the 
target (so that the net effect of the evolution is to replace a single-pomeron 
exchange, cf. Fig. 2. a, by a double-pomeron exchange, Fig. 2.d). ii) In Fig. 
2.f, the same vertex describes gluon merging in the target wavefunction. 

By 'iterating' the diagrams in Fig. 2, it is straightforward to deduce the di- 
agrammatic interpretation of the equation satisfied by (T^)y for any N > 1. 
Although this seems like a trivial exercise, we nevertheless exhibit the cor- 
responding results for (T^)y in Figs. 4.b and c (projectile evolution) and 
4.d and e (target evolution), only to emphasize that these are not all the di- 
agrams which would be expected within perturbative QCD to the order of 
interest. (Fig. 4. a shows the corresponding tree-level diagram.) This lack of 
perturbative completeness is to be related to our previous remarks concerning 
the lack of saturation effects in the projectile evolution and, respectively, the 
lack of gluon number fluctuations in the target evolution. As we shall shortly 
see, both omissions refer in fact to the same physical process, only seen from 
different points of view. 

Specifically, in Figs. 4.f and g we show two diagrams which are formally 
of the same order in a s as the diagrams for linear evolution in Figs. 4.b and 
d, yet they were not included in the previous analysis. For the diagram in 
Fig. 4.f, this absence is easy to explain: this diagram is subleading at large 
N c , so it is naturally missed by Eq. (3.1); as a matter of facts, this diagram is 
correctly included in the full version of Balitsky equations (valid for arbitrary 

7 Indeed, according to Eq. (3.6), the four gluons which enter this vertex from the 
above are pairwise coupled into two color singlets, so like the two gluons emerging 
from the vertex; after subsequent evolution, any of these 2-gluon singlet exchanges 
would be eventually converted into BFKL ladders, or 'pomerons'. 
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(a) 



(b) 




C < 




e ? 



(d) 




(e) 





(g) 

Fig. 4. Diagrams for dipole pair scattering: (a) tree—level contribution; (b, c) one 
step projectile evolution; (d, e) one step target evolution; (f) a diagram suppressed 
at large N c ; (g) the missing diagram, which is of leading order in both ot s and N c . 

N c ), and also in the JIMWLK formalism. As for the second diagram, Fig. 4.g, 
this is easily recognized to describe the incriminated process, which is missing 
in either Balitsky, or JMWLK, picture of the non-linear evolution: Seen from 
the perspective of projectile evolution, this diagram describes the interaction 
between two dipoles in the projectile, that is, a saturation effect, whereas from 
the perspective of the target, it rather corresponds to gluon splitting. 



By inspection of the diagrams in Fig. 4, one can understand what was 
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the argument for ignoring such contributions in the JIMWLK evolution of the 
target, and also when is this argument expected to break down. Consider the 
dependence of the processes in Fig. 4 upon the gluon density in the target: A 
2-point function counts like the density, a 4-point function, like the density 
squared, etc. Then the standard, BFKL, diagram in Fig. 4.d is proportional 
to the gluon density squared, the recombination diagram in Fig. 4.e scales like 
the third power of the density times an additional factor a s , while the splitting 
diagram, Fig. 4.g, scales likes a s times the density 8 . Clearly, this last contri- 
bution is indeed negligible whenever the density is large enough. For instance, 
in the saturation region, the density is 0(1/ a s ), so both the linear diagram, 
Fig. 4.d, and the recombination diagram, Fig. 4.e, contribute on equal foot- 
ing, whereas the splitting diagram, Fig. 4.g, is totally irrelevant. However, the 
situation changes in the dilute regime at high transverse momenta, where the 
density is very low (the average gluon occupation number can be much smaller 
than one); then, all diagrams are suppressed except for the splitting one, Fig. 
4.g, which shows how a 4-point correlation gets generated via the splitting 
of a 2-point one. Clearly, this is the dominant mechanism for constructing 
higher-point correlation functions in the dilute regime. 

To summarize this discussion, the Balitsky- JIMWLK equations include 
a single type of triple pomeron vertex, which corresponds to either pomeron 
splitting, or pomeron merging, depending upon the perspective that we have 
over the evolution. Still, both merging and splitting are necessary in order 
to have a complete description of the approach towards saturation, which 
is consistent with perturbative QCD. Among the approaches to high energy 
evolution proposed so far, the dipole picture includes only splitting (and thus 
can be used to describe particle number fluctuations at high-/cj_), whereas 
the JIMWLK evolution contains only merging (so it describes correctly the 
saturation effects at high density). As explained in the Introduction, our goal is 
to combine the respective virtues of these two approaches in order to achieve 
a more complete description of the high-energy evolution of the scattering 
amplitudes. Our strategy can be summarized as follows: 

Since the physical regimes in which the interesting phenomena — sat- 
uration and fluctuations (or merging and splitting) — play important roles 
are kinematically well separated, it should be possible to use a 'piecewise' 
description in which each regime is covered by a different formalism. Such a 
description should correctly capture the essential physics that we are inter- 
ested in, although it would probably fail to reproduce the fine details of the 
transition between the two regimes. Specifically, in what follows we shall stick 



Overall, there is the same power of a s associated with the diagrams in Figs. 
4.d and 4.g. But when interpreting these diagrams in terms of gluon densities, one 
should remember that one factor of a s is included in the definition of the gluon 
density. 
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to target evolution, and we shall rely on the JIMWLK formalism in the high- 
density regime, and on the dipole picture in the dilute one (which means, in 
particular, that we shall limit ourselves to the large- N c approximation). Thus, 
both gluon merging, Fig. 4.e, and gluon splitting, Fig. 4.g, will be encoded in 
the target wavefunction, although their respective theoretical descriptions will 
be different. But, as we shall see, these descriptions can be naturally merged 
with each other in the equations for the scattering amplitudes, thus leading to 
a generalization of the Balitsky-JIMWLK hierarchy which includes the effects 
of gluon number fluctuations in the dilute regime, and thus generates pomeron 
loops 9 through iterations. 



4 A toy model borrowed from statistical physics 

A main difficulty in our present analysis is the lack of a unified theoretical 
framework for describing both saturation and particle number fluctuations 
in the hadron lightcone wavefunction. It is therefore instructive, before we 
embark ourselves in a more detailed study of QCD, to make a detour through 
a simple particle model borrowed from statistical physics which contains both 
physical ingredients (particle number fluctuations and recombination) in a 
unified setting. This model is formulated similarly to the QCD dipole picture 
that we shall discuss in the next section — there is a system of particles 
whose distribution evolves in time according to a master equation — but it 
goes beyond the dipole picture by including recombination effects (so, in that 
respect, it may be viewed as a non-linear version of the dipole picture). 

For this model, we shall study the effects of fluctuations on the evolution 
equations for particle number correlations, and, in particular, we shall show 
that it is possible to exactly reformulate the dynamics as a specific Langevin 
equation. This is interesting since the Langevin equation is better suited for 
numerical studies than the original master equation. Inspired by this, we shall 
later be able to also reformulate the equations for scattering amplitudes in 
QCD as a corresponding Langevin equation. Besides, some of the conceptual 
and technical issues that we shall clarify in this section will appear again later, 
in the context of QCD. 

For more clarity, in this section we shall consider the simplest version 
of the statistical model: a zero-dimensional system (no spatial dimension) 
in which the only interesting distribution is that of the particle number. In 

9 The simplest 'pomeron loop' is a one-loop diagram generated by a (BFKL) 
pomeron which splits into two which then recombine with each other. Thus, for 
this diagram to exists, one needs the simultaneous presence of triple pomeron ver- 
tices for splitting and merging. 
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the Appendix, we shall then discuss a more refined model (taken from Refs. 
[38,43]), which has one spatial dimension, and which in spite of its relative 
simplicity leads to a Langevin equation very similar to the one that we shall 
find in QCD. This is the sFKPP equation [39,40] (and Refs. therein). 

Specifically, in this section we consider a system whose state is described 
by the number n of particles A. One can think of all these particles as being 
located at some unique lattice site, in which case n represents the occupation 
number of that site. The dynamics is determined by two random elementary 
processes: under an infinitesimal step dY in time, a particle can split into two 
with a rate a : 

A^>A + A, (4.1) 
while two particles can recombine into one with a rate 2/3 : 



A + A^A, 



(4.2) 



Clearly, the probability P(n) = P(n, Y) to find a configuration with n particles 
changes with time, and the corresponding master equation reads 



dP(n) dP{n) 



+ 



dP{n) 



dY 



dY dY 

where the first term in the r.h.s. is the contribution of splitting: 
dP(n) 



(4.3) 



dY 



= a[(n- l)P(n - 1) - nP(n)] , 



(4.4) 



while the second one is due to recombination 
dP{n 



dY 



(3 [n{n + l)P{n + 1) - n(n - l)P{n)\ . 



(4.5) 



The expectation value of an observable 0(n) is given by 

(0(n)) Y = J2P(n,Y)0(n), (4.6) 

n 

and its evolution equation can be obtained in general from Eq. (4.3)-(4.5). 
For example, it is not hard to find how the average values of the "number 
operators" (n), (n 2 ),... change with time. We have 



d{n) 

~w 

d(n 2 ) 
dY 



a(n) — (3 (n 2 ) — (n) , 

2(n 2 ) + (n)] -p\(n 3 }-3(n 2 } + (n) 



■ a 



(4.7) 
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We shall also need the normal-ordered number operators, defined as 

:n K := n(n - 1) • • • (n - k + 1). (4.8) 

To understand why we call this a "normal-ordered product" , imagine a second- 
quantized formalism in which particles of type A are created, or annihilated, 
by standard Fock space operators, a) and, respectively, a, with [a\ a] = 1. 
Then the particle number operator is n = aJa, and we have n 2 = o)ao)a = 
o)a)aa + a) a =:n 2 : +n, so :n 2 := n{n — 1), so like in Eq. (4.8). Note further- 
more that : n 2 : gives zero when acting on a state with a single particle (unlike 
n 2 , which yields one). Thus, in the present context, the normal-ordering in 
: n 2 : is introduced in order for this operator to properly count the pairs of 
particles (so that, e.g., a single particle cannot be counted like a pair), and 
similarly for :n K : . Moreover, as we shall see in the next section, and also in 
the Appendix, the use of normal ordering is the only way to avoid equal-point 
singularities for systems with spatial dimensions. 

The evolution equation for = (:n K :) takes a rather simple form: 



dn^ 
~~dY 



= a 



Kn w + k(k 



l)n 



(«-i) 



P 



Kn (K+l) + k(k 



l)n 



(«) 



(4.9) 



But the simplicity is rather illusory, as is coupled in the evolution with 
n( K_1 ) and rS K+1 \ so Eq. (4.9) is just a particular equation within an infinite 
hierarchy. In what follows, we shall study various aspects of this hierarchy: In 
particular, we shall discuss the respective roles of the multiplication and re- 
combination processes, and also the relative importance of fluctuations versus 
mean field aspects at different stages of the evolution. But before doing that, 
let us first introduce a different stochastic process which is mathematically 
equivalent to the one that we have studied so far, in the sense of generating 
exactly the same correlations as obtained by solving the hierarchy in Eq. (4.9). 



Specifically, let us consider the following Langevin equation 

dbTb i 

— = an- ph 2 + y/2 [ah- f3h 2 } v = A + Bv , 



(4.10) 



where v(Y) is a Gaussian white noise : (u(Y)) = and iy{Y)v(Y')) = 5(Y — 
Y f ), and the right hand side must be understood with the Ito prescription 
for time discretization (recall that a Langevin process is not differentiable): 
Namely, if one writes Y = je, where j — 0, 1, 2, . . . and e is the length of the 
time step, then Eq. (4.10) should be properly understood as 



n 



71,- 



a 7i,- 



Pn 2 + X l2 



a 



(3ti 



Aj + Bj v 



(4.11) 



with 
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te> = 0, (u t u J ) = -J lJ , (4.12) 

which shows that updating the variables n at time-step j + 1 only requires 
knowledge of the variables at time-step j. The noise term in Eq. (4.10) is said 
to be multiplicative since it is multiplied by a function of n. It is easy to see 
that h = is an unstable fixed point for the evolution described by Eq. (4.10), 
whereas h = a/ {3 is a stable fixed point, and corresponds to saturation. If 
the initial condition satisfies < n(0) < a/ (3, this will remain true for all 
subsequent times. 

For an arbitrary function F(n), Eq. (4.10) implies the following evolution 
equation 

d -^^ = (AF'(n)) + 1 -(B^F"(n)}. (4.13) 

By taking F(n) = n K , it is straightforward to show that the hierarchy emerging 
from the Langevin problem is indeed equivalent to the hierarchy of the toy 
model, Eq. (4.9), provided we identify (n K ) with n^ K \ In particular, the noise 
term in Eq. (4.11) is responsible for the fluctuation terms in Eq. (4.9), i.e., the 
terms proportional to k(k — 1) within the square brackets. 

It is in fact easy to check directly at the level of the original hierarchy, 
Eq. (4.9), that = (a/ (3) K is indeed a fixed point of the evolution. Moreover, 
in the presence of fluctuations, we do not expect other fixed points (see also 
below). Thus, for sufficiently large Y, the occupation number will saturate 
to n sat = a /(5 independently of the initial conditions. We would like this 
saturation value to be large, n sat ^> 1, in order for the analogy with QCD to 
be as close as possible. We shall therefore choose a»^, which implies that 
the recombination effects (proportional to (3) will be parametrically suppressed 
as compared to the growth effects (proportional to a) except in the saturation 
regime — a situation similar to QCD. Thus, if one starts in a dilute regime 
{n <C n sat ) at Y — 0, then it is possible to follow the evolution in the early 
stages (prior to saturation) by using the linearized version of Eq. (4.9), as 
obtained after neglecting recombination : 



dY 



= an 



(i) 



dY 



= 2a 



(4.14) 



As we shall see in Sect. 5, the above equations are the analog of the equa- 
tions satisfied by the dipole densities in the dipole picture. Note that these 
equations are still coupled with each other, because of the fluctuation terms: 
The term linear in in the r.h.s. of the evolution equation for describes 
a fluctuation in which one particle splits into two (so this is the analog of 
the dipole splitting term in Eq. (2.2)). One may expect fluctuations not to be 
important at large times where ^> but in general this is not right. 
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To see this, consider the solution to Eqs. (4.14) with generic initial conditions 

n ( 1 )(0) = n and n( 2 )(0) = n£ ] : 



n {1) {Y) =n Q exp(aY), 



n (2) (Y) = [n{, 2) + 2n J exp(2aY) - 2n exp(aY), . . . (4.15) 

An interesting initial condition is no = 1 and = for n > 2, meaning that 
there is a single particle in the initial state. (In QCD, this would correspond to 
the case where one starts the evolution with a single dipole, or, more generally, 
with an isolated fluctuation at high k±.) Then Eqs. (4.15) show that so long 
as aY < 1, dominates over in the r.h.s. of the second equation (4.14), 
and thus is driving force for the growth in the number of pairs. In the absence 
of fluctuations, the pair number would remain zero for ever. For larger Y, 
starts to dominate over (as it increases faster), but the ratio 

n {2) (Y) 

' " ' -2-2exp(-aY), (4.16) 



[nW(Y)]< 

approaches 2, and thus differs from the naive prediction n^^Y) xs [n^(Y)] 2 
of the mean field approximation. In fact, one can check on Eqs. (4.15) that 
the only way to satisfy this mean-field factorization at intermediate values 
of Y (prior to saturation) is to assume a large occupation number already 
in the initial conditions: no ^> 1, n 2 ^ ~ njj ^> no, etc. (In QCD, this would 
correspond, e.g., to starting the evolution with a very large nucleus at Y — 0.) 

This conclusion is not altered by the addition of the recombination terms, 
which merely provide saturation: For the model at hand and for generic initial 
conditions (no ~ 0(1)), the mean field approximation is justified only in the 
saturation regime achieved at very large Y : Y ^ Y ~ (1/a) ln(a//5). For 
a more general system which involves also spatial dimensions, so like QCD 
or the reaction-diffusion model to be studied in the Appendix, particles can 
escape from the bulk, so there will be always a region in phase-space in which 
the density is low and the mean field approximation is not applicable, however 
large is Y and independently of the initial conditions (even if one starts with 
a large nucleus): this is the tail of the distribution, in which the evolution is 
driven by fluctuations. 

In view of the comparison with QCD, it is useful to notice that the analog 
of the (dipolar) Balitsky-JIMWLK equations in the present toy model are the 
equations obtained after neglecting the fluctuation terms in Eq. (4.9), that is 



K 



an {K) - (3n {K+1) (no fluctuation) . (4.17) 



dY 

(Alternatively, the analog of Eq. (3.1) - - the generating equation of the 
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Balitsky-JIMWLK hierarchy — is the deterministic equation obtained after 
ignoring the noise term in the Langevin equation (4.10).) Neglecting fluctua- 
tions in the recombination term is indeed harmless, since this term is important 
only in the high-density regime, where fluctuations become irrelevant anyway. 
But for the growth term, which acts also in the dilute regime, the fluctuations 
are crucial, as shown by the previous analysis. In particular, the reduced hier- 
archy in Eq. (4.17) could never generate a high density system from an initial 
state which involves only one particle. 

Another effect of the fluctuations is to wash out some remarkable prop- 
erties of Balitsky equations, like the existence of exact factorized solutions 
[41,42] : Eq. (4.17) can be solved with the Ansatz = c K [n^] K provided 
the coefficients c K satisfy c K = c£ _1 , in agreement with Ref. [42]. (In particu- 
lar, = (l/c 2 )(a/(3) K is a one-parameter family of fixed points for Balitsky 
equations.) But it can be easily checked that such a factorized solution does 
not exist for the complete equations including fluctuations, that is, Eq. (4.9). 

More generally, we expect the fluctuations to wash out any sensitivity 
to the initial conditions after a sufficiently large evolution. For the zero- 
dimensional model in this section, this is simply the statement that at large Y 
the solution will converge to an unique fixed point, namely = (a/(3) K , in- 
dependently of the initial state. For more complicated systems, which involve 
also spatial dimensions, we expect universality at late times (in the sense of 
insensitivity to the initial conditions) also for more complex aspects of the 
dynamics, like the behavior of high energy scattering amplitudes in QCD. 



5 Fluctuations and evolution in the dipole picture 

We now return to QCD with a discussion of particle number fluctuations 
in the dipole picture. To that aim, we shall not rely on the original formula- 
tion of the dipole picture due to Mueller [5,6], but rather on its alternative 
formulation in Ref. [32] (see also [37]), which is better suited for a study 
of fluctuations. The main difference between these two formulations is that, 
loosely speaking, 'they put the evolution at different ends' (see also Figs. 5 
and 6 below). More precisely, in Mueller's original formulation, the evolution 
proceeds via splitting at the highest rapidity end: the rapidity increment dY 
is used to accelerate the original dipole (the one which has initiated the evolu- 
tion), which then undergoes an additional splitting, whose effects propagate in 
the whole configuration. But in this picture the correlation between the high- 
rapidity dipole which has split and the low-rapidity one that we measure is 
distributed over many steps of rapidity, and thus is difficult to trace back. 
By contrast, in the formulation in Ref. [32], the splitting occurs in one of the 
dipoles at the lowest rapidity end, which was itself generated in the previous 
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step of the evolution, and which is the direct parent of the dipole we measure. 
This makes it easier to follow correlations associated with splitting, which are 
particularly important in the first few steps after the splitting occurs. 

Following Ref. [32] , the system of dipoles generated by the evolution up to 
rapidity Y of an original dipole with coordinates x and jjo will be described 
as a stochastic ensemble of dipole configurations endowed with a probabil- 
ity law which evolves with Y according to a master equation. Specifically, a 
given configuration is specified by the number of dipoles N and by N — 1 
transverse coordinates {z^} = {zi, z 2 , ...Zn-i}, such that the coordinates of 
the N dipoles are (z ,Zi), (z±, z 2 ),...,(zn-i, Zn), with z = x and z N = y . 
The probability P^{{zi}\ Y) to find a given configuration obeys the following 
evolution equation (this is similar to Eq. (4.4)) : 

dP N (z 1 ,...z N _ 1 ;Y) a s \^ 

*N{Zl, ■■■Zn-1, J- ) 



dY 2tt 



J "d 2 ''z M(zi- U Zi, z) 



i=i 

+ 77^ Z_] M(Zi-i, Zi+i, Zi) Pjv_i(Zi, Zi-i, Z i+ i, Zat„ : 

271 i=l 

(5.1) 

where we have defined the shorthand notation for the dipole kernel 

M(x,y,z)= {X ~ y)2 (5.2) 
(x — z) 2 (y — z) 2 

The expectation value of an operator O which depends only the dipole sizes 
is given by 

oo „ 

(0(Y)) = ^ dr N P N ({z i };Y)0 N ({z i }), (5.3) 

N=l J 

where the phase space integration is simply air N = d 2 Zid 2 z 2 . . . d 2 ZN-i- Then 
by using the master equation (5.1) one can show that 



d(0{Y)) _ a_^j dfN Pn({Zi} . y) £ J d 2 zM(Zi _ uZuZ) 

[-0 N ({ Zi }) + N+l ({ Zi ,z})], (5.4) 



dY 

ui ^„ N=1 . 



X 

where the z argument in O^+i is to be placed between and Zj. 

In what follows, we shall use Eq. (5.4) to derive evolution equations for the 
dipole number densities. Because of the difference alluded to before with re- 
spect to the original formulation of the dipole picture, the final equations that 
we shall obtain are not the same as the corresponding equations in the early 
literature on the dipole picture (see Refs. [5,6,30,31]), although the two sets of 
equations are mathematically equivalent (in the sense of providing identical 
results for identical initial conditions). Very recently, Levin and Lublinsky [8] 
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have obtained the same equations that we shall derive here, but they used a 
more straightforward approach, based on a generating functional [37], which 
however hides some of the subtle points that we shall emphasize below, and 
which in our opinion are important to properly understand the results. 

Consider first the average dipole number density. The corresponding op- 
erator for a iV-dipole configuration is 

N 

n N (x,y) = £<5 (2) (^-i - x)S^( Zj - y). (5.5) 
i=i 

With On = un, the expression within the square bracket in the second line 
of Eq. (5.4) becomes 

A i (x 1 y,z) = -S^(z i . 1 -x)S^(z i -y) 

+ <5 (2) (zi-i - x) 5^{z -y) + 6<®(z - x) 5^{z t - y). (5.6) 

After simple manipulations we arrive at the following evolution equation for 
the average of the dipole number density n Y (x,y) = (n(x,y)) Y ■ 



dn Y (x,y) _a L , 2 
<)Y 2~ ' " 



M(x,y,z)n Y (x,y) 

M(x, z, y) n Y (x, z) + M(z, y, x) n Y (z, y) 



= ^ j d 2 zlC X y Z ®n Y (x,y). (5.7) 

As anticipated, this is not the same as the standard equation for the dipole 
number density in the dipole picture. In the latter, the splitting occurs in 
the original dipole (x ,y ), so the more complete notation n Y (x, y \x , y ) is 
needed for the dipole density. Then the standard equation reads 



dn Y (x,y\x ,y ) a 



dY 



ol r r 
= 2tt J d2zM ( x 0'y0' z )[ ~ n Y (x,y\x ,y ) 



+ n Y (x, y\x , z) + n Y (x, y\z, y ) , (5.8) 



and is also recognized as the dipole version of the BFKL equation [1]. The 
two equations (5.7) and (5.8) are illustrated in Figs. 5 and 6, respectively. In 
spite of the formal differences, these equations are nevertheless equivalent 10 , 
as shown in Ref. [32]. 

Now we turn to the dipole pair density. As in the toy model considered in 
the previous section, the corresponding operator for a given iV-dipole config- 



10 This can be checked by using the symmetry property (x — y) 4 n Y (x, y\xQ,y ) 
(x - yo) 4 n Y (x ,y \x,y). 



27 




Fig. 5. The one-step evolution of the average dipole number density as described 
byEq. (5.7). 




Fig. 6. The same as in Fig. 5, but from the perspective of Eq. (5.8). 
uration is defined as the normal-ordered product of single densities: 

(2) 

n y N > (xx, yi, x 2 , y 2 ) = : n N (x 1 , y x ) n N (x 2 , y 2 ) ■ (5.9) 

= n N (x 1 ,y 1 )n N (x 2 ,y 2 )-5 {2) (x 1 - a; 2 )5 (2) (?/i - y2)n N (x 1 ,y 1 ). 

In fact, it is precisely this normal-ordered operator that is automatically re- 
produced by the generating functional used by Mueller [5,6], and also by Levin 
and Lublinsky [8,37], but the importance of normal ordering has not been em- 
phasized before, nor the difficulties arising when using the simple product (as 
opposed to the normal-ordered one) have been discussed. 

The use of normal ordering is required by both physical and technical 
considerations 11 . On the physical side, note that Eq. (5.9) is equivalent to a 
restricted double sum: 

N 

n%\x uyi] x 2l y 2 ) = £ gVfa-! - xx)8 { - 2 \z j - y 1 )8 { - 2 \z k _ 1 - x 2 )8^(z k - y 2 ), 

j,k=l 
3+k 

(5.10) 

which does not allow for "pairs" made of the same dipole. On the technical 
side, note that the subtraction performed in the second line of Eq. (5.9) elim- 
inates the 'ultraviolet' singularity of the simple product n^(x\,y\)n^(x 2 ,y 2 ) 
in the limit where the two dipoles have identical coordinates 12 . Because of 

11 We are grateful to Larry McLerran for his insightful observations on the role of 
normal ordering in this particular context. 

12 This is manifest on the expression in Eq. (5.10) which is well defined when x\ = x 2 
and y\ = y 2 . The only dangerous term in this limit would have been the term with 
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that, the evolution equation for the normal-ordered pair density is well defined 
in the continuum limit, unlike for the simple product. 



Specifically, by using Eq. (5.9), we see that the square bracket in the second 
line of Eq. (5.4) becomes (recall also Eq. (5.6)) 



n 



(2) 
JV+l 



n 



(2) 
N 



Aj(xi, z) n N (x 2 , 2/2) + Aj(a3 2 , 2/2, z)n N (x 1 , 2/1) 



(5.11) 



+ Aj(cci, 2/1, z) Ai(x 2 , 2/2, z) - 5 {2) (x 1 - x 2 ) 5 (2) (yi - 1/2) \(x u 2/1, z). 



The contribution of the first two terms to the evolution equation is rather easy 
to obtain and reads (to simplify writing, we use the notation introduced in 
Eq. (5.7) together with n^' = (n^) 



xiyiz 



n ( $\x 1 ,y 1 ;x 2 ,y2) + linear, 



(5.12) 



where the linear terms not explicitly shown arise from normal-ordering the 
product tin(xi, 2/1) n^{x 2 ^ y 2 ). By themselves, these linear terms are divergent 
when X\ = x 2 and/or 2/1 = y 2 , but in the complete calculation they exactly 
cancel the singular terms coming from Aj(cci, 2/1, z) Ai(x 2 , 2/2, z) and from the 
last term in Eq. (5.11). All the remaining terms are non-singular, and add up 

(2) (2) 

to give the following contribution to n N+1 — n N : 

'5^(z i _ 1 -x 1 )5^(z-y 1 )5^(z-x 2 )5^(z i -y 2 )} + {l ^ 2}, (5.13) 

which in turn gives the following contribution, linear in the dipole number 
density n Y , to the evolution equation for the dipole pair density : 



^-M(x 1 ,y 2 ,x 2 )n Y (x 1 ,y 2 )5 (2 \x 2 -y 1 ) + {l <-> 2}. 
Putting everything together we finally arrive at 
dn { y\x 1 ,y 1 ]x 2 ,y 2 ) a s 



dY 



2^ 



/ 



d 2 z!C 



n { Y ] {xi,yi]X 2 ,y 2 ) 



xiyiz 



+ M(x u 2/2, x 2 ) n Y {x u 2/2) 5 {2 \x 2 - 2/1) 



(5.14) 

(5.15) 
{1^2}. 



(2) 

The terms proportional to n Y in the r.h.s. of the above equation describe 
the normal BFKL evolution of the pair density, while the terms linear in n Y 
are recognized as fluctuations in which the two measured dipoles arise from 
the splitting of the same parent dipole (and thus are contiguous with each 
other). These are precisely the fluctuations discussed in Sect. 2. As explained 
there, and also for the statistical model in Sect. 4, the main effect of these 
fluctuations is to drive the growth of pair density in the dilute regime where 
n Y dominates over n Y ^ in the r.h.s. of Eq. (5.15). 



j = k in the unrestricted double sum. 
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Fig. 7. The one-step evolution of the dipole pair density n Y as described by 
Eq. (5.15). The virtual BFKL term is not shown. 

The evolution described by Eq. (5.15) is illustrated in Fig. 7. Once again, 
this equation is not the same as the standard equation for n Y which appears 
in the literature [5,6,30,31] (but it coincides with the corresponding equation 
of Levin and Lublinsky [8]). The standard equation is rather obtained by 
putting the evolution in the original dipole (so like in Eq. (5.8)), and its r.h.s. 

(2) 

involves n Y together with terms quadratic in ny, but not also terms linear in 
ny. Fluctuations are therefore more difficult to recognize on that equation. 

For later use it is also convenient to express the fluctuation term in terms 
of the dipole sizes = x, L — yi and the corresponding impact parameters 
bi = (xi + i/i)/2, instead of the coordinates of their quark and anti-quark 
legs. One thus finds 

a s (n + r 2 ) 2 / r 2 \ (2) / r l + r 2 \ r 

(5.16) 

which is the same as Eq. (2.2), as expected. 

Note finally that via the same manipulations as above it is straightforward 

(k) 

to construct the equation obeyed by the /c-body dipole density n Y within the 
dipole picture. This is the analog of the a-piece of Eq. (4.9) in the statistical 

(k) 

model, that is, it involves a piece proportional to n Y which represents the 
normal BFKL evolution, and a piece linear in fiy -1 ' which describes fluctua- 
tions and dominates at low density. The corresponding equations can be found 
in Ref. [8]. 



6 Fluctuations and saturation in dipole— CGC scattering 

In the previous section, we have seen that the dipole picture provides a 
natural theoretical framework to describe gluon number fluctuations in the 
dilute regime, and for large N c . In what follows, we shall exploit the relation 
between dipole densities and scattering amplitudes in order to transfer the 
effects of fluctuations from the target wavefunction to the evolution equations 
for the amplitudes. In the dilute regime, where the fluctuations are important, 
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the scattering amplitudes are simply proportional to the dipole densities in 
the target (see Sect. 6.1), and it is therefore straightforward to translate the 
previous equations in Sect. 5 into corresponding equations for the amplitudes. 
The ensuing equations can be then extended to the high-density regime by 
simply adding the same non-linear terms as in the Balitsky-JIMWLK equa- 
tions (since the fluctuations cannot change significantly the non-linear aspects 
of the dynamics). These manipulations will lead us, in Sect. 6.2, to a new set 
of equations which include both fluctuations and saturation, and which are 
the main result in this paper. Finally, in Sect. 6.3, we shall show that, mediat- 
ing a coarse-graining in impact parameter space, this infinite set of equations 
can be replaced by a single Langevin equation, which looks formally like the 
Kovchegov equation supplemented by a noise term. 

6.1 From dipole densities to scattering amplitudes 

Let us start with the scattering between a single external dipole ("the 
projectile") and a dilute target whose wavefunction is described in the dipole 
picture. The very statement that "the target is dilute" does not refer to the 
target alone, but rather is a statement about how the target looks like on 
the resolution scale of the projectile, which is its transverse size. Namely, this 
implies that the size r of the external dipole is small as compared to the local 
saturation length l/Q s (Y,b) in the target at the dipole impact parameter 
b. Under these conditions, the scattering with the external dipole acts as 
a measure of the density of gluons (or dipoles) with transverse size r in the 
target wavefunction. This is so since the dipole-dipole scattering is quasi-local 
in transverse phase-space, as we recall now: 

A priori, the external dipole can scatter off any of the dipoles present in 
the target wavefunction, so its scattering amplitude can be evaluated as : 

T(r, b,Y)=J d 2 n J d 2 b x T (r, n, b - b{) n(n, b u Y), (6.1) 

where n(ri,bi,Y) is the number density of dipoles with size V\ and impact 
parameter b\ in the target wavefunction, as produced after a rapidity evolu- 
tion Y, and T (r, r±, b — bi) is the amplitude for the scattering between two 
elementary dipoles. Eq. (6.1) holds for a given configuration of the target wave- 
function, that is, for a given scattering event. The corresponding expression for 
the average amplitude is obtained after averaging over all the configurations 
in the target, as explained in Sect. 5. 

Note that in writing Eq. (6.1) we have restricted ourselves to a single- 
scattering approximation, as appropriate in the dilute regime. This equation 



31 



is correct only so long as T(r,b,Y) <C 1. For consistency, it is sufficient to 
evaluate T to lowest order in perturbation theory, that is, to 0(a 2 ), which 
corresponds to the exchange of two gluons. The corresponding expression is 
well known in the literature, but here we shall need only its general properties. 
Namely, % ~ a 2 when the two dipoles have similar sizes and are relatively 
close to each other in impact parameter space, but it decreases very fast, as 
7o(t*,t*i,6 — bi) ~ a 2 r 2 r\/{b — foi) 4 , for a large separation \b — b\\ ^> r> 
between the two dipoles (with r> = max(r, ri)). One can thus replace bi — > b 
as the argument of the density n in the r.h.s. of Eq. (6.1), and at the same 
time evaluate the integral over 6 X as: 

d 2 b x T (r, n, 6 - In) ~ y r> rf 2 6 r ~ a*r* , (6.2) 

where b r = b — bi and r< = min(r, 77). 

To evaluate also the remaining integral over r 1; one needs to know the 
corresponding dependence of the dipole density. So long as we are in the 
genuine BFKL regime, characterized by a non-trivial "anomalous dimension" , 
the dominant behavior of n(ri, b, Y) is expected to be of the form 

n(r u b,Y) ~ ±(rlQl(b,Y)y°, (6.3) 

where 70 ~ 0.63, 1 — 70 is the anomalous dimension, and Q 2 (b, Y) oc e XoisY is a 
line of constant density (up to the kinematical factor 1/rf) in the ln(l/ri) — Y 
plane. (This line is parallel to the saturation line.) Strictly speaking, Eq. (6.3) 
holds for the average density, but for sufficiently large Y, we expect the dipole 
distribution to self-average even in a single event, and therefore show the same 
gross features as the average distribution. After inserting Eq. (6.3) in the r.h.s. 
of Eq. (6.1), one can easily check that the ensuing integral is saturated by 
T\ ~ r, so, by dimensional arguments, the final result reads simply: 

T{r,b,Y) - a 2 s r 4 n(r,b,Y) . (6.4) 

In the dilute BFKL regime in which we are interested here, this equation is 
exact up to some numerical fudge factor of order one, that we do not control 
in the above approximation, and which should depend upon the value of the 
anomalous dimension. 

Eq. (6.4) can be immediately extended to the scattering of two or more 
external dipoles. In a particular event, the scattering amplitude for a system 
of several external dipoles is simply the product of the individual amplitudes 
for each dipole. For instance, 
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T Y 2) (n, b i; r 2 , b 2 ) = T(n, h, Y) T(r 2 , b 2 , Y) 

~u\r\r\ :n(r 1 ,b 1 ,Y)n(r 2 ,b 2 ,Y): , 



(6.5) 



where the final expression is correct up to some unknown fudge factor. Note 
that this final expression involves the normal-ordered dipole pair density, as 
defined in Eqs. (5.9)-(5.10). This is so since we neglect here the possibility that 
both external dipoles scatter off the same dipole in the target. As explained 
in Sect. 2, this is a good approximation so long as the external dipoles have 
different sizes and/or impact parameters. 

6.2 The equations for the scattering amplitudes 

By using Eqs. (6.4) and (6.5) together with the equations for the average 
one-body and two-body dipole densities in the target, as established in Sect. 
5, one can immediately deduce the equations satisfied by the average scattering 
amplitudes within the present approximations. 

For the scattering amplitude of a single dipole {T(r,b)) Y = (T(x,y)) Y , 
we thus find the BFKL equation, as expected: 



Note that after multiplication by the fourth power of the dipole size, the 
various dipole kernels appearing in the original equation (5.7) for the dipole 
density have been all converted into the same kernel A4(x, y, z). 

Eq. (6.6) applies so long as T <C 1 . But we know already what should be 
the corresponding generalization to the high density regime where T ~ 1 : 
In that regime, the target wavefunction is described by the CGC formalism, 
which should not be significantly altered by fluctuations, because the latter are 
relatively unimportant when the density is large. Thus, the general equation 
for (T) can be simply obtained by adding to Eq. (6.6) the appropriate non- 
linear term generated by the JIMWLK evolution of the target. As discussed 
in Sect. 3, this non-linear term is (— T(x, z)T(z, y))y- With this addition, 
Eq. (6.6) becomes identical to the first equation (3.3) in the Balitsky hierarchy. 

Consider now the scattering amplitude for two dipoles (T^ (xi, y±; x 2 , y 2 ))y 
= (T(xi,yi)T(x 2 ,y 2 ))Y- In the dilute regime, this is obtained by combining 
Eqs. (6.5) and (5.15), while in the high density regime this should also include 
the non-linear terms (involving (T^) Y ) induced by the JIMWLK evolution. 



d(T(x,y)) Y 
dY 



2^ 



I 



d 2 zM(x, y,z)(- T(x, y) + T(x, z) + T(z, y)) y 




(6.6) 
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By putting all these ingredients together, we find the following evolution equa- 
tion (with the concise notation introduced in the second line of Eq. (6.6)) 



Among the various terms in the r.h.s. of Eq. (6.7), those in the first two lines 
were already present in the corresponding Balitsky equation, while the term 
linear in (T)y in the third line (that we have written in the mixed notations 
of Eq. (2.2)) is a new term, which takes into account the effect of fluctuations. 
The unknown fudge factor k ~ 0(1) in this term has been introduced to 
parametrize our uncertainty concerning the precise relation between T and n, 
cf. Eqs. (6.4) and (6.5). Note that this uncertainty does not affect the other, 
more standard, terms in the r.h.s. of Eq. (6.7), which come out the same as 
in the Balitsky-JIMWLK equations. 

Let us discuss the properties of this new term in some detail. First, this 
is of order a 2 s (T), so it is parametrically suppressed when T a 2 (since in 
that case already the uncorrelated piece (T)(T) of (T^) is much larger than 
the fluctuation term). In particular, the new term is certainly unimportant in 
the saturation regime, where Eq. (6.7) reduces to the corresponding Balitsky 
equation. On the other hand, the fluctuation term dominates in the very dilute 
regime where (T) <C a 2 , and thus is responsible for the growth of (T^) in 
that regime. For instance, if we start with a target made of a single dipole 
at Y = 0, then (T^) = 0, and the rise of (T^) in the early stages of the 
evolution is driven by the fluctuation term. A similar discussion applies to a 
rare high-A;^ (or small size) fluctuation generated in the wavefunction of an 
arbitrary target. We conclude that the usual BFKL evolution is not applicable 
in the dilute regime, in contrast to what one could naively expect. This failure 
will be further analyzed in Sect. 7. 

Consider also the geometry of the fluctuation term in the transverse plane, 
as manifest on Eq. (6.7). As anticipated in Sect. 2, for the external dipoles to 
feel the effect of fluctuations, they must be contiguous with each other, i.e., 
x 2 = Vi or 2/2 = x i- Moreover, for contiguous dipoles of unequal transverse 
sizes, the importance of the fluctuation is controlled by the size of the smallest 
dipole: r\r\j(r\ + r 2 ) 2 ~ r\. This is so because of the geometry of dipole 
splitting (see Fig. 1) together with the fact that, in order to scatter, two 
dipoles have to overlap with each other. When a (target) dipole splits into 



d(T^(x l ,y 1 ;x 2 ,y 2 )) Y 
dY 
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d 2 z (M x1v1z <g> (T {2 \x 1 ,y 1 -x 2 ,y 2 )) Y 
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Fig. 8. Diagrammatic illustration of the fluctuation term in Eq. (6.7) : the original 
dipole with size r\ + splits at the time of the interaction into two new dipoles 
with sizes r\ and, respectively, r2, which then scatter off two external dipoles. 




Fig. 9. Two steps in the evolution of the average scattering amplitude of a single 
dipole: the original amplitude (left) and its evolution after two steps (right). 

two dipoles of very different sizes, the small child dipole is necessarily located 
near the edge of the parent dipole. Thus, for the scattering to take place, 
the impact parameter of the incoming small dipole should be located within 
a distance r< from the edge of the parent dipole (of size |t*i + r 2 \ ~ r>). 
This condition introduces the geometrical penalty factor r</r> manifest in 
the fluctuation term in Eq. (6.7). To conclude, the effects of the fluctuations 
are larger for incoming dipoles of comparable sizes. 

Consider also the diagrammatic interpretation of the fluctuation term. 
This is illustrated in Fig. 8, where we represent the target as a collection of 
dipoles (as appropriate in the dilute regime), one out of which splits into two 
dipoles which then interact both with the external ones. This diagram is the 
dipole picture analog of the 2-to-4 splitting diagram in Fig. 4.g. Of course, 
the dipole picture cannot be used to also describe the non-linear effects in 
the previous equations, which correspond to saturation effects in the target 
wavefunction. To interpret the latter, we should rather resort to the CGC 
picture. In Fig. 9, we depict a particular two-step evolution of the single 
dipole amplitude (T), as emerging from Eqs. (3.3) and (6.7). The upper part 
of this diagram is the same as the splitting process in Fig. 8. The lower part 
describes the recombination of the four gluons resulting from that splitting 
into two; this is a diagram of the CGC formalism (compare to Fig. 2.f) which 
is encoded in the non-linear term in Eq. (3.3). The overall diagram in Fig. 9 
represents the simplest pomeron loop generated by the present equations. 
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It is straightforward to write down the generalization of Eq. (6.7) valid 
for an amplitude (T^) with k > 3. This will be the corresponding (dipolar) 
Balitsky equation supplemented by fluctuation terms linear in (T^ k ~^) which 
account for the possibility that a pair of dipoles (r^rj) within (T^) which 
are contiguous with each other be generated through the splitting of a unique 
dipole of size r*j + Vj from (T^ -1 )). The ensuing hierarchy of equations, which 
generalizes the Balitsky hierarchy at large N c by including the effects of gluon 
number fluctuations in the dilute regime, represents our main result in this 
paper. As already manifest on Eq. (6.7), the equations in this hierarchy are 
afflicted with an uncertainty concerning the normalization of the fluctuation 
terms, which reflects the approximation (6.4) used in their derivation. As we 
shall explain in Sect. 7, this uncertainty does not affect the leading order pre- 
dictions of these equations at very large Y and small a s , but it may influence 
the subleading effects. 



6.3 The Langevin equation 



The equations for scattering amplitudes constructed in the previous sub- 
section are similar to, although more complicated than, the general equations 
(4.9) for particle number correlations in the statistical toy model. In Sect. 4, 
we have seen that the hierarchy in Eq. (4.9) can be equivalently represented 
by a Langevin equation with an appropriate noise term. In what follows we 
shall show that, mediating some approximations which allow one to get rid 
of the impact parameter dependence of the amplitudes, a similar Langevin 
equation can be written also in QCD. 

The Langevin equation will be seen to emerge naturally once that we suc- 
ceed to rewrite the equations for the scattering amplitudes in a form quasi-local 
in b. To that aim, we need to assume that the target is quasi-homogeneous 
over distances of the order of the size of the external dipoles. This is a reason- 
able approximation so long as we are interested in the local evolution towards 
sat ur at ion/unit arity, i.e., in the local 'blackening' of the target. In particular, 
this is sufficient for a study of the energy dependence of the saturation mo- 
mentum. On the other hand, this approximation cannot be used to study the 
transverse expansion of the black disk, or to compute total cross-sections. 

Under this assumption, all the terms in the previous equations except for 
the fluctuation terms become local in b. For instance, the r.h.s. of Eq. (3.3) 
for (T(x,y)) = (T(r,b)) involves T(x,z) = T(r',b+^) (with r' = x-z), 
which is approximated as T(r',b). Then, Eq. (3.3) becomes: 
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d(T(r,b)) Y _ a s r r 2 

dY ~ 2irJ r«(r-r') 2 1 J 

- T(r, 5) + T(r', b) + T(r - r', b) - T(r', b)T(r - r', b)) y , 

which is indeed local in b. Similarly, all the standard terms in the equation 
for (T^fVi, 61; r 2 , b 2 )) — that is, the terms in the first two lines in the r.h.s. 
of Eq. (6.7) — become bilocal in 61 and b 2 . But the fluctuation term there 
remains non-local, because of the 5-functions imposing contiguity. We have 
indeed (compare to Eq. (5.16)) 



9(T( 2 )(r 1 ,6 1 ;r 2 ,6 2 



Y 



dY 



= «<^ ( r i ra \ 2 (nri+r 2 ,bi); y 

fluct . 2tt (n + r 2 ) 2 \ iy 

(bi-bi- 7 ^) + {1-2}, (6.9) 



where we have approximated T(t*i + r 2 , 61 — w T(ri + r 2 , 61). 



At this level it becomes natural to perform a Fourier transform to mo- 
mentum space, by introducing 



= I^T2 e ~ lkrT (r,b). (6.10) 

Up to an overall normalization of 0{\/6t s ), <p{k, b) can be interpreted as the 
gluon phase-space occupation number, or the 'unintegrated gluon distribution' 
[26,44]. The benefit of using the momentum space is twofold: (i) After Fourier 
transform, the non-linear terms in the evolution equations (e.g., the term 
quadratic in T in Eq. (6.8)) become local in k. (ii) Since the fluctuation term 
(6.9) depends upon the transverse sizes r*i and r 2 only through their sum 
T\ + 7*2, its Fourier transform will be diagonal in momentum. 

Specifically, Eq. (6.8) yields in momentum space : 



(6.11) 

where the linear part is recognized as the BFKL equation for the unintegrated 
gluon distribution, while the non-linear term provides saturation. 

Furthermore, the Fourier transform of the fluctuating contribution to 
(T( 2 )), Eq. (6.9), is obtained as 
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d (ip(k 1 ,b 1 )ip(k 2 ,b 2 ))^ 



dY 



= a s 2Ka 2 s 5 {2 \k l -k 2 ) (6.12) 

fluct. 



where the factor of 2 appears because the two terms in the r.h.s. of Eq. (6.9) 
have given identical contributions. Note that the exponential in the last inte- 
gral selects r < 1/fci, thus the fluctuation term is important only for transverse 
separations \b\ — b 2 \ ^ 1/fci between the two dipoles. 

To also rewrite the fluctuation terms in a form which is local in b, we shall 
proceed to a coarse-graining. This is physically motivated, since an external 
dipole of size r cannot distinguish details in the target on a transverse scale 
much smaller than r. That is, the gluon (or dipole) distribution measured 
by the external dipole is necessarily averaged in b over a disk of radius r 
around its impact parameter. After Fourier transform to momentum space, 
this corresponds, by the uncertainty principle, to a coarse-graining over an 
area ~ 1/k 2 . To implement this, we shall divide the impact parameter space 
into cells of area 1/k 2 , and average the scattering amplitudes over any such a 
cell. Consider a particular cell with center at bj. Averaged quantities within 
that cell are then defined as: 



ip^k) = k 2 [ d 2 bip(k,b), ip 2 (k) = k 2 I d 2 bcp 2 (k,b), (6.13) 

etc., where Sj(fc) ~ 1/k 2 is the area of the cell. In fact, given our previous 
assumption that the target is quasi-homogeneous over the area covered by the 
dipole, it is clear that the averaging above is tantamount to simply replacing 
(p(k, b) everywhere within a cell by its value at the center of that cell: (pi(k) ~ 
tp{k,bi). Therefore, (fi 2 (k) [(pi(k)] 2 , so, by itself, this coarse-graining does 
not introduce additional correlations. 

The coarse-graining is therefore trivial for the standard (BFKL + non- 
linear) terms in the evolution equations: it amounts to replacing tp(k, b) — > 
(Pi(k) (for a dipole hitting the target in cell i) in any of those terms. The 
non-locality of the BFKL kernel poses no special difficulties for this coarse- 
graining since, e.g., the integral over p in Eq. (6.11) is dominated by momenta 
p ~ k. Then Eq. (6.11) implies that the equation for the average amplitude 
(tpi(k)) in cell % involves the double-scattering amplitude within the same cell, 
(ip 2 (k)). To compute this, it is therefore enough to study the scattering of two 
external dipoles which hit the target within the same cell. This is also the 
interesting case for a study of fluctuations: as discussed after Eq. (6.12), the 
fluctuations vanish when the external dipoles fall in different cells. 

Consider therefore the equation satisfied by ((pi(ki)ipi(k 2 )) . The coarse- 
graining is non-trivial only for the fluctuation term, in which case it implies 
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(note that ki = k 2 in this term, cf. Eq. (6.12)): 



-•tL/^f^^WrM)) -*?«*<*)>• (6.14) 

The integration over &2 has been used to remove the 5-function (whose support 
is indeed included in £j(/ci)), while the integration over bi has served together 
with Eq. (6.10) to construct (ipi(k)). 

Let us summarize here the first two equations of the hierarchy, as applying 
to the scattering within cell i: 



(6.15) 

and, respectively, 



%(%(fc2))y - \ [d 2 p k\ I ( p 2 \ \ 

-(^ 2 (fcl)^(fc2)) y , 

+ / t a s 2 ^(fc 1 -fc 2 )^(^(fc 1 ))y} 

+ {l^2}. (6.16) 

Clearly, these equations are (quasi) local in b in the sense that they involve 
only operators pertaining to cell i. Needless to say, this property holds also for 
the higher equations in the hierarchy. Within the present, quasi-homogeneous, 
approximation, different cells in the impact parameter space evolve indepen- 
dently from each other. 

It is now straightforward to identify the Langevin equation which gener- 
ates the same correlations as the hierarchy at cell %. This reads (from now on, 
we omit the cell index i, as there is no coupling between different cells): 



+ ^2Ka 2 s ip(k)u(k) (6.17) 
where v(k,Y) is a Gaussian white noise : (u(k,Y)) = and 
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Mfci^iM*^)) = -5(Y 1 -Y 2 )5^(k 1 -k 2 )k 2 l . 



(6.18) 



Eq. (6.17) must be understood with a rapidity discretization prescription of 
the Ito type, as explicitly shown in Eq. (4.11). 

Eq. (6.17) should be compared to Eq. (3.1) which, we recall, is the gen- 
erating equation of the Balitsky hierarchy. Clearly, Eq. (3.1) is the same as 
the (Fourier transform of the) deterministic part of Eq. (6.17), but the lat- 
ter involves also the noise term responsible for fluctuations. In that respect, 
Eq. (6.17) is closer to the Langevin equation (4.10) for the statistical toy 
model. Note, however, that in contrast to Eq. (4.10), the noise term in the 
corresponding equation (6.17) for QCD does not vanish at saturation. This 
reflects the fact that, in our previous analysis of QCD, we did not include 
the effect of fluctuations on the recombination terms (an effect which would 
go beyond the dipole picture that we have used to describe fluctuations in 
QCD). However, this omission is truly harmless: In the high density regime 
at/near saturation, fluctuations are anyway unimportant, and there is no need 
to suppress them explicitly 13 . In fact, it is rather obvious on Eq. (6.17) that 
the effects of the noise term are important only so long as ip(k) ^ a 2 . 



7 Physical discussion 



The Langevin equation (6.17) turns out to be the natural generalization 
to QCD of the stochastic Fisher-Kolmogorov-Petrovsky-Piscounov (sFKPP) 
equation (see the Appendix), to which it reduces in a standard approxima- 
tion to the BFKL kernel known as the 'diffusion approximation'. This will be 
discussed in Sect. 7.1. Then, in Sects. 7.2 and 7.3, we shall use some known 
results about the sFKPP equation [38,45,46] (see also Refs. [39,40] for recent 
reviews and more references) in order to explore the physical consequences of 
Eq. (6.17) for QCD. In particular, in Sect. 7.2, we shall make the connection 
with the physical discussion and the results in Ref. [34]. Then, in Sect. 7.3, 
we shall discover a rather dramatic consequence of the fluctuations, namely 
the breakdown of the BFKL approximation in the regime where the average 
gluon density is small. 



13 Numerical simulations within the context of statistical physics [38,45] have con- 
firmed that the sFKPP equation, in which the noise term has as a coefficient 
\J aip(l — ip) — with a< 1 and tp = 1 at saturation — and the so-called Reggeon 
model, where the coefficient of the noise is simply y/ct<p, lead indeed to identical 
results for the measured correlations. 
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7.1 Relation with the sFKPP equation 



In order to understand the relation between our Langevin equation Eq. (6.17) 
in QCD and the sFKPP equation in statistical physics, it is useful to first re- 
call a few facts about the deterministic part of Eq. (6.17), as obtained after 
neglecting the noise term there. If we also neglect correlations in the initial 
conditions, then this deterministic part is simply the momentum-space version 
of the BK equation, which can be more compactly rewritten as: 



dy(p(p,Y) = a s x(-d p )(f -a a <p (7.1) 

where p = lnk 2 /kl (with k some arbitrary momentum scale of reference), 
xd) — 2-0(1) — 0(7) — "0(1 — t) i s the Mellin transform of the BFKL kernel, and 
x(— dp) is an integro-differential operator defined via the series expansion of 
x(l) ( see below). In a series of papers [29], Munier and Peschanski have argued 
that Eq. (7.1) is in the same universality class as the Fisher-Kolmogorov- 
Petrovsky-Piscounov (FKPP) equation, which appears as a mean field ap- 
proximation to a variety of stochastic problems in chemistry, physics, and 
biology [39,47]. To understand this correspondence, notice that in the di- 
lute regime at high transverse momenta, k 2 3> Q 2 S (Y) or p > p s (Y), with 
p s (Y) = hi Ql(Y)/k%, we have <C 1, so the dominant behavior of ip(p,Y) is 
determined by the linear part of Eq. (7.1), which is the BFKL equation [1]. 
One thus finds that the dominant dependencies upon p and Y can be isolated 
out into an exponential factor: 



(p(p,Y) = e-^ p ~ XoiisY) ^(p,Y), (7.2) 
where 70 and Aq are pure numbers determined by the BFKL kernel as [2] : 



7ox'(7o) = x(7o), A = ^, (7.3) 

7o 

which implies 70 ~ 0.63 and Ao = x'ilo) ~ 4.88, and the function ip(p, Y) is 
comparatively slowly varying, so its behavior can be studied by using a limited 
expansion of the operator x{~ dp) around x(lo) '■ 



Xi-dp) = x(7o) + x'(lo)(-d p - 70) + lx"(lo)(-d p - 7o ) 2 . . . 

^-X d p + D (-d p - l0 ) 2 , (7.4) 

where D = x"(7o)/2 and in writing the second line we have also used 
Eq. (7.3). The approximation which consists in keeping only the terms to 
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second order in this expansion, as explicitly shown in the second line above, is 
generally referred to as the "diffusion approximation" , and is equivalent to a 
saddle point approximation to the solution to BFKL equation in Mellin space, 
v In the case of BFKL equation, one can use Eqs. (7.2) and (7.4) to show that 
ip(p, Y) obeys a diffusion equation. For the non-linear BK equation (7.1), one 
can argue that the main effect of the non-linear term is to introduce an ab- 
sorptive boundary condition on the diffusion equation for ip [27]. Alternatively, 
one can study the full non-linear equation which emerges within the diffusion 
approximation: 



dyip(p, Y) = a s ( - \ d p + D (-d p - 7o) 2 )y? - a s y? 2 . (7.5) 

Munier and Peschanski have observed that, up to a linear change of variables 
and an appropriate rescaling of </?, Eq. (7.5) is the same as the FKPP equation: 

d t u(x, t) = dlu(x, t) + u(x, t) (f - u(x, *)) , (7.6) 

which at large times describes an uniformly translating front which propagates 
from the stable state u — 1 into the unstable one u — and decays exponen- 
tially at large x (far ahead the front) [39,47]. In QCD, this front is already 
visible on Eq. (7.2), which at large Y describes a traveling wave located at 
p = p s (Y) ~ \oa s Y which has an exponential slope 70 and propagates in 
'time' a s Y with uniform velocity Ao- (This behavior is the origin of geometric 
scaling [35] for the BK equation [26,27].) By using known properties of the 
FKPP equation [39] or, alternatively, by solving the diffusion equation for ip 
with an absorptive boundary condition [27], one can deduce the dominant cor- 
rections to the velocity and the shape of the front due to non-linearities. For 
large Y and p 3> p s (Y), one finds ip(p, Y) ~ (p — p s (Y))/Y 3 ^ 2 , and therefore: 

dpsiY) 3 f 

Returning to the original Langevin equation (6.17) and assuming that the 
dominant p-behavior of the solution is still given by the exponential e^ op 
(which turns out to be right indeed), one can use again the diffusion approxi- 
mation (7.4) to deduce a simplified form of the equation (r = a s Y): 



d T (p(p,r) = -Xod p <p + D (-d p - 7 ) 2 (p - </ + ^2k a 2 s ip is(p, r), (7.8) 
with the Gaussian white noise : 



(u(p, r)> = 0, (u(p, T )v(p>, t')> =-S(t- t') S(p - p') . (7.9) 

71 
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Up to a simple change of variables, the equation above is essentially the same 
as the sFKPP equation [40]: 



d t u(x,t) = d 2 x u + u(l-u) + ^—u(l-u)u(x,t). (7.10) 

In this analogy, the density of particles per site u(x,t) in the stochastic par- 
ticle model corresponds to the scattering amplitude (or unintegrated gluon 
distribution) ip(p, r) in QCD, and 1/a 2 plays the same role as N (the number 
of particles per site at saturation), which is reasonable since we recall that 
n ~ l/cKg is also the dipole occupation number at saturation. Therefore, the 
mean field approximation emerges in the limit where the occupation numbers 
are large, i.e., N — > oo for the particle model and a 2 — > in QCD. 

Note finally the limitations of this correspondence, as inherent in the use 
of the diffusion approximation. First, Eq. (7.8) predicts that behind the front 
(f saturates at a constant value (p = 7q-Do (or m = 1 for Eq. (7.10)), whereas 
the original Langevin rather yields ip(k,Y) as ha(Q s (Y)/k) [i.e., ip(p,Y) m 
(p s (Y) — p)/2] for p <C Ps{Y), as most easily seen by using Eq. (6.10) to- 
gether with the fact that T(r) = 1 for r ^> 1/Q S (Y). Second, for fixed Y and 
sufficiently large p (essentially, such that p — p s (Y) > p s (Y) [26]), the BFKL 
'anomalous dimension' 1 — 7 approaches to zero, which signals the transition 
to a regime dominated by the DGLAP dynamics. Clearly, in this regime it is 
not possible to expand the BFKL kernel around the saturation exponent 70, 
as we did in Eq. (7.4). Thus, we do not expect the simplified equation (7.8) 
to describe correctly the transition to the DGLAP regime. 

7.2 Some results from sFKPP equation and their consequences for QCD 

Even if somehow simpler than the original Lagevin equation (6.17), the 
sFKPP equation (7.8) (or (7.10)) remains complicated, because of the simul- 
taneous presence of the non-linear term for recombination and of the multi- 
plicative noise term. Fortunately, this equation has been extensively studied 
in relation with problems in statistical physics — mostly through numerical 
simulations, but also via some analytical methods — , with results that we 
shall briefly describe here and then adapt to the QCD problem at hand. In 
applying these results to QCD, one should however keep in mind the possible 
limitations of the correspondence between Eqs. (6.17) and (7.10), as mentioned 
at the end of the previous subsection. 

It has been rigorously demonstrated [46] that the front generated by the 
sFKPP equation (7.10) is compact : For any t, there exists a x r (t) such that 
u(x, t) = for x > x r (t), and also a Xi(t) such that u(x, t) = 1 for x < 
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xi(t). One expects a similar property for Eq. (6.17), although in that case 
the transition to the saturation region behind the front may not be as sharp, 
because the coefficient of the corresponding noise term does not vanish at 
saturation. 



Numerical simulations show that for large but finite N, the front generated 
by Eq. (7.10) propagates with an asymptotic velocity which is smaller than 
the corresponding velocity v q = 2 for the FKPP equation (7.6). Moreover, with 
increasing N the convergence of vn towards vq is extremely slow: vq — vn ~ 
1/ In 2 N when N >> 1. For the corresponding QCD problem, this implies the 
following asymptotic F-dependence of the saturation scale : 

Xs = Km d fM. ~ A - - 2 .f. when a s « 1, (7.11) 

dr In (1/a 2 ) 

(recall that r = a s Y). Thus, with decreasing a s , the saturation exponent X s 
converges only slowly towards the respective mean-field value in Eq. (7.3). 

Since the noise term in Eq. (6.17) is important only for relatively small 
ip a 2 , we expect the shape of an individual front (in its comoving frame) 
to be rather well described by the solution (7.2) to the deterministic equation 
(7.1) at all the points where p 3> a 2 (or u 3> 1/N for Eq. (7.10)). This is seen 
indeed in the numerical simulations. One can then estimate the width of the 
front as follows: It is the noise term in Eq. (6.17) which abruptly cuts down the 
growth of tp ahead of the front, but this requires tp to be as small as a 2 . Since in 
the front region behind the tip ip behaves like ip ~ e -70 ^ - ^, c f jt^ (7.2), we 
conclude that ip decreases from 1 to a 2 over a range p — p s ~ (I/70) ^(l/ 1 ^ 2 )) 
which should be a good estimate for the width of the front. 

However, because of the fluctuations inherent in the noise term, different 
realizations of the same evolution will lead to an ensemble of fronts which all 
have the same shape, but are displaced with respect to each other along the 
p-axis. That is, the position p s of the front is itself a random variable, charac- 
terized by an expectation value (p s (r)), which for large r increases according 
to Eq. (7.11) (since this is the common asymptotic behavior of all the fronts 
in the ensemble), and also by a dispersion a 2 = (p 2 ) — (p s ) 2 , which is expected 
to rise linearly with r : (T 2 (t) ~ D {t t (since the front executes a random walk 
around its average position). The numerical simulations to the sFKPP equa- 
tion confirm this behavior, and show that the front diffusion coefficient Df r 
scales like 1/ln 3 iV when N ^> 1. For QCD, this in turn implies: 



V 

3/777^7 when ( 7 - 12 ) 



m d (l/« 2 ) 

which vanishes, as expected, when a s — > 0, but only very slowly. 



44 



Note that the results (7.11) and (7.12) are only logarithmically sensitive 
to the coefficient of 0(a 2 ) of the noise term in Eq. (6.17), or (7.8). Thus, 
the leading order estimates in the limit a s — > 0, as shown in Eqs. (7.11) and 
(7.12), are not affected by our uncertainty concerning the fudge factor k in the 
noise. On the other hand, the next-to-leading order correction already will 
be sensitive to the precise coefficient under the log, and thus to k (since, e.g., 
l/ln 2 (a//«* 2 ) w 1/ ln 2 (l/a 2 ) - 2 ln(a//c)/ m 3 (l/a 2 )). 

The results of the sFKPP equation alluded to above turn out to be con- 
sistent with numerical and analytic studies of the stochastic particle models 
whose mean field approximation is the FKPP equation (7.6). This agreement 
should not come as a surprise: As shown in the Appendix, the sFKPP equa- 
tion emerges precisely as the continuum description of such discrete particle 
models. In particular, in that context, Brunet and Derrida [48] have given a 
simple heuristic argument which explains the 1/ In 2 N scaling of the velocity 
correction vq — at large N, and also allows one to compute the corre- 
sponding coefficient C. Namely, they have observed that, in the presence of 
discreteness, diffusion should replace local growth as the main mechanism for 
front propagation: Indeed, in order for the growth term (the linear term in 
u(x, t) in Eq. (7.6)) to be effective, there must be at least one particle per bin 
(or lattice site). Thus, the only way that a particle can move to an originally 
empty bin ahead of the front is via diffusion from the previously occupied bins 
on its left. To mimic that, Brunet and Derrida proposed a modified determin- 
istic equation obtained by inserting a cutoff 6{u — l/N) in the growth term in 
Eq. (7.6) (since a density u ~ 1/N corresponds to a site occupation number of 
(9(1)). A simple analysis of this equation then implies [48] vq — v^ — C/ In 2 N 
for N >> 1, with a value for C which is indeed consistent with the numerical 
studies of both particle models and the sFKPP equation. After translation to 
the QCD problem of interest here, this in turn implies: 



A s ^ A - '7 ' when a s « 1. (7.13) 

2 In (l/a 2 s ) 

It so happens that the QCD coefficient C = 7r 2 7 x"(7o)/2 is numerically large, 
C 150, so the corrective term in the equation above can be trusted only for 
extremely small values of a s , which are physically unrealistic. To our knowl- 
edge, there is no analytic argument allowing one to understand the scaling 
(7.12) of the front diffusion coefficient, or to compute the coefficient T> there. 
(See however the discussion in Sect. 4 of Ref. [40].) 

Furthermore, in the context of QCD, the above results (7.11) and (7.12) of 
the sFKPP equation corroborate the conclusions obtained in Ref. [34] through 
an analogy between the high-energy problem in QCD and some specific parti- 
cle models in statistical physics. This in turn demonstrates that the Langevin 
equation (6.17), or, more generally, the stochastic equations presented in Sect. 
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6.2, provide the correct evolution law underlying the physical picture put for- 
ward in Ref. [34]. By (numerically) solving these equations, one can now go 
beyond the analysis in Ref. [34] and study the evolution for realistic (non- 
asymptotic) values of Y and a 2 s . 

7. 3 Front diffusion and the breakdown of the BFKL approximation 

It has been argued in Ref. [34] that, as a consequence of the diffusive 
wandering of the front, the 'geometric scaling' property characteristic of the 
individual fronts — i.e., the fact that a particular front realization propagates 
as a travelling wave, T(p, Y) m T(p — p s (Y)) (cf. Eq. (7.2)), so that its shape 
in the comoving frame is not changed under the evolution — is actually bro- 
ken after averaging over the statistical ensemble of fronts generated by the 
stochastic evolution up to large Y. In what follows, we shall demonstrate that 
the consequences of the front diffusion are in fact even more dramatic, as they 
entail the breakdown of the BFKL approximation in the dilute regime, where 
this approximation is usually assumed to work 14 . What we mean by that more 
precisely is that, at very high energies, the standard BFKL equation (i.e., the 
linearized part of Eqs. (3.3) or (7.1)) fails to correctly describe the evolution 
of the average scattering amplitude (T(p)) Y even in the regime where this 
amplitude is small, {T(p)) Y <C 1. This is so because, in the presence of fluc- 
tuations and for sufficiently large Y, average quantities like (T) Y or (T^) Y 
are dominated by those fronts within the statistical ensemble which are at 
saturation for the values of p of interest, and this even when p is well above 
the average saturation momentum {p s ) Y , where (T(p)) Y <C 1 indeed. 

Note that in this subsection we return to the notation T (rather than 
ip) for the amplitude, that is, we prefer to work in coordinate space, where 
p = ln(l/r 2 /co). Then, for a given front realization, with saturation momentum 
p s (Y), the scattering amplitude saturates to T = 1 behind the front, and takes 
the scaling form T(p, Y) ~ e -7o(p-p s ) w itfiin a finite range in z ahead of the 
front. (Recall that the front is compact.) For the present purposes, it suffices 
to describe this behavior with the simple interpolation: 



which neglects the compact nature of the front; this is harmless since, as we 
shall see, the tail of the distribution at large z does not contribute to average 
quantities in the regime of interest. 

Because of the diffusive nature of the front wandering, the values of p s are 
14 We are grateful to Al Mueller for helping us clarifying this point. 




for p < p s 
for p> p s 



(7.14) 
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distributed according to the probability density 



P(Ps) = 



exp 



(Ps - (Ps)Y 



(7.15) 



where (p s ) ~ X s a s Y and a 2 ~ D{ r a s Y, as discussed before. Then the average 
amplitude (T) is determined by 



(T(P, (Ps))) 



dp s P(Ps)T(p,p s ), 



(7.16) 



and higher-point correlations can be computed similarly. For simplicity, in 
what follows we shall restrict ourselves to higher correlations evaluated at 
equal points, e.g., (T^ 2 \p,p)) = (T 2 (p)). Then, one can obtain (T n ) by simply 
replacing 7o — > n 7o in the subsequent formulae for (T). 



By using Eqs.(7.14) and Eq.(7.15) it is straightforward to show that 



(T) = iErfcg) + iexp(^- 7o ; 



2-Erfc(*-^ 



(7.17) 



where z = p — (p s ) and Erfc(x) is the complimentary error function, for which 
we recall that 



Eric(x) 



2 - 



exp(— x 2 ) 
y/nx 



exp(— x 2 ) 
^/nx 



for x <C — 1 
for x = 

for x ^> 1. 



(7.18) 



For what follows it is important to notice that the two terms in the r.h.s. of 
Eq.(7.17) arise from the saturating piece and, respectively, the exponentially 
decaying piece, of Eq.(7.14). 

The behavior of (T) as a function of z depends upon the competition 
between a (the width of the Gaussian distribution of the fronts) and I/70, 
which characterizes the exponential decay of the individual fronts. Since a 
grows like \/Y whereas 70 ~ 0(1), one may conclude that the typical situation 
at high energy is such that a 3> l/ 7 o- But the diffusion coefficient (7.12) 
vanishes when a s — > 0, even though only slowly, so for sufficiently small a s 
at fixed Y one can imagine also the situation where a <^ l/ 7 o- So, we shall 
consider both cases here: 

• a < I/70 

In this regime, the front diffusion should play no role, and this is indeed 
what we find. Specifically, when z <^ —a, by making use of Eq.(7.18) one sees 
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that the first term in Eq.(7.17) gives 1, while the second one is negligible. On 
the contrary, when z ^> a, the first term is small while the second one gives 
exp(— Jqz). Thus, not surprisingly, one finds that the average amplitude (T) 
retains the single event profile, except in the short interval \z\ < a where it gets 
smoothed. The BFKL dynamics, characterized by the anomalous dimension 
7o, is still visible, and the mean field approximation holds since (T n ) ~ {T) n . 

• a > I/70 

In this regime, the correlations are dominated by fluctuations, and the 
BFKL behavior is washed out everywhere except at extremely large distances 
ahead of the front (where however the present approximations cannot be 
trusted since, e.g., the compact support property of the front, and also higher 
correlations in the diffusion of the front, may play a role). 

Specifically, for all values of z such that z <C 70O" 2 one finds that (T) is 
dominated by the first term in Eq.(7.17): 

(T) ~ l - Erfc (^j for - 00 < z « 7o cr 2 , (7.19) 

which, as mentioned earlier, is the contribution from the saturating pieces 
of the single events, and thus is independent of 70. This estimate holds, in 
particular, in the range a <C z <C 7 <t 2 where (T) is small, (T) <C 1, yet very 
different from the corresponding BFKL prediction. In particular, there is no 
trace of geometric scaling, in agreement with Refs. [33,34]. 

In fact, within the whole range in z in which Eq.(7.19) is valid, the higher 
correlations (T n ) are given by this same expression, that is 

(T n ) ~ (T) for - 00 < z < 7o cT 2 , (7.20) 

which signals a total breakdown of the mean field approximation, except in 
the saturation regime where (T) ~ 1. One can even find a window within 
which (T n ) varies very slowly, namely: 

(T n ) ~ (T) ~ ^ for \z\ < a. (7.21) 

The different behaviors encountered when increasing a (and thus Y) are 
illustrated on the example of (T) and (T 2 ) in Fig. 10. 

Let us conclude this section, and also the paper, with a final remark con- 
cerning the role of fluctuations in the high-energy evolution in QCD, and, in 
particular, on their interplay with saturation. It has been a recurrent theme in 
this paper that fluctuations in the dilute regime act as a seed for the growth 
of higher-point correlations which then, through their subsequent evolution, 
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Fig. 10. Evolution of (T) and (T 2 ) with increasing a. 

generate the non-linear terms responsible for saturation. This is the (more 
or less) expected part of the scenario, in which fluctuations and saturation 
are important at different ends of the spectrum, and the amplitude T(p, Y) is 
reasonably well described by the mean field approximation (the BK equation) 
everywhere except at the tip of the distribution where T is very small: T < a 2 . 
This picture is indeed what emerges from the dynamical equation (6.17), but 
only so long as one considers a single front realization (i.e., a single event). If, 
on the other hand, one considers the statistical ensemble of fronts — which is 
what one needs to do in order to compute average quantities — , then the dis- 
persion of the fronts due to fluctuations has a rather unexpected consequence, 
which is to render the correlation functions sensitive to saturation even in the 
(formally) weak scattering regime, where the average scattering amplitude is 
small, and thus wash out the normal BFKL behavior. 
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A A one— dimensional reaction-diffusion model 

In this Appendix we shall study a much more sophisticated model, than 
the one we did in Sect. 4, which shares a lot of similar features with QCD. 
We consider a system whose state is described by the number of particles A 
at every site % in a one-dimensional lattice; each state is of the form {n^} = 
{ni,ri2, ■■■,rii, ...}, where rii is the number of particles at site i. The dynamics 
contains three possible processes: 

• A particle can split locally into two at a rate a, i.e. 

A^Ai + Ai. (A.l) 

• Two particles can recombine locally into one at a rate 2/3 , i.e. 

A. + A.^A,. (A.2) 

• A particle can diffuse to a neighboring site at a rate /i, i.e. 

Aj — > Aj + i, Aj — > Aj_i. (A.3) 

Our task is to derive the evolution equation for (rij) and (riiiij) from the 
master equation, which we do not write explicitly here, but the corresponding 
changes in probabilities will be obvious from the analysis below. Then we will 
study the continuous limit, which will lead us to the necessity of requiring 
normal-ordering, and we will end this Appendix by introducing the, relevant 
to the problem, Langevin equation. 

Since splitting and recombination take place locally, their contribution to 
the time evolution of (rij) is easy to obtain. For both these processes, one can 
have either gain or loss under a step dY, depending on the number of parti- 
cles in the initial configuration. There are two terms arising from the splitting 
part dP({rii}; Y)/dY\ a (all the other splittings which occur at lattice points 
different than the one we measure, i.e. the i-th one, cancel each other): 
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o Loss — arti P(..., rii, Y) 

o Gain +a (rii — 1) P(..., rii — 1, Y), 



and two terms from the recombination part dP({rii}; Y)/dY\ /3o 

o Loss -Porii(ni- l)P(...,n h ...;Y) 

o Gain +[3 (m + 1) m P(...,rn + 1, Y). 



To calculate d(rii)/dY \ a ,/3 , we appropriately shift the n's wherever needed to 
obtain 



d{rii) 



dY 



(A.4) 



= J2P(M;Y)F aA) ({n t }), 

a,(3 {m} 

where, after very simple algebra, one can find that F a ^ ({rii}) = ani—j3 (n% — 
rij). Therefore we arrive at 



d(rii) 



dY 



a(rii) - (3 



(A.5) 



Now let us also study in detail the evolution of (rii) due to diffusion. Under 
a step dY in time, the diffusion can drive an particle from site i to % ± 1. 
Again, in all cases we can have either gain or loss depending on the number 
of particles in the initial configuration. Overall there are eight terms arising 
from the diffusive part dP({ni}; Y)/dY\ fJi : 



o Loss 



- urii P(..., rij_i, rii, n i+1 , Y) 

- iin i P(...,n i - 1 ,n i ,n i+1 , Y) 

- Iirii-! P(..., rii_i, rii, n i+1 , Y) 

- nn i+ i P(...,n i - 1 , rii, n i+ i, Y) 



from % — > % — 1 

from i — > i + 1 

from % — 1 — > i 

from i + 1 — > i 



Gain 



+ /i (nj + 1) P(...,nj_i - + l,n i+1 , Y) from i — > i — 1 

+ ii (rii + 1) P(...,ni-i,rii + - 1, Y) from i — > i + 1 

+ /i(n i _i + l)P(...,n i _i + l,n i -l,n i+ i,...;y) from i-l^i 

+ /i + 1) P(..., rij-i, nj — 1, n i+ i + 1,...;Y) from i + 1 — > i 



As before, in order to calculate d{rii) /dY]^, we appropriately shift the n's 
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wherever needed to obtain 

d(ui) 



dY 



= ^P(K};F)F,(K}), 

At {rii} 



(A.6) 



with F M ({nj}) = n{n i+ i + nj_i — 2nj). Therefore we arrive at 

d(rii) 



dY 



(A.7) 



where, for reasons that will become obvious shortly, we introduced the short- 
hand notation 

V 2 rii = n i+ i + ni-x - 2n i+1 . (A. 8) 

Now, when considering the pair density we need to study the three separate 
cases (nf), {riiUi+i) and (niiij) for > 2. Here we shall skip the derivation, 
which can be done by following the same steps of the (n») case, and only give 
the final set of the evolution equations which read 



diyiiUj) 



dY 

d{nirij) 



dY 

dfaiiij) 



= a(2n i n j + 5^), 



(A.9) 
(A.IO) 



dY 



jjL(.niV 2 nj + niV 2 rij 



+ %(V 2 ni + Arii) - {Sij-i + 6ij +1 )(rii + 



n. 



(A.11) 



Let's turn our attention to the continuous limit. Let the lattice spacing be 
equal to A, the coordinate x = iA and the particle single and pair densities 
n(x) = rii/A and n(xi)n(x2) = fiifij/A 2 . In order to be economical in the 
notation we shall drop the bars immediately. Then the single density equation 
becomes 



d{n(x)) 

dY =a(n(x))-p A 



(n 2 {x)) - l(n(x)> 



+ /.A 2 + 0(vl A 4 ). 

(A.12) 



Notice that the next to last term justifies the shorthand notation adopted 
in Eq.(A.8). Considering first this diffusion term, we naturally impose that 
/j,A 2 = D — fixed in the A — > limit. Then all the higher order terms can be 
ignored. Similarly, we require f3 A = (3 — fixed for the recombination term. It 
is only the recombination process which generates potential singularities in the 
continuous limit, since it involves the pair density. However, even though the 
second term in the square bracket becomes divergent, we expect the first term 
to do so in the same limit in such a way to cancel the singularity. Indeed, the 
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whole bracket corresponds to the X\ 
pair density 



x 2 limit of the average normal-ordered 

n (2) (xi,x 2 ) = (:n(xi)n(x 2 ) :) = (n(xi)n(x 2 )) - - x 2 )(n(xi)), (A.13) 



and thus we finally arrive at 

dn<V(x) 



dY 



an 



x) - f3n [2 \x,x) + £> 



d 2 n^(x) 
dx 2 



(A.14) 



with the obvious identification n^(x) = (n(x)). Of course, in the absence of 
recombination, these subtleties would have not appeared. But they do appear 
for all the processes, i.e. even for splitting and diffusion, when we try to take 
the continuous limit in the equation of the pair density. As an illustration, let's 
concentrate on the diffusion contribution. When we consider the \i — j\ > 2 
case we arrive at an equation with no singularities, since all the terms in the 
second line of Eq.(A.ll) vanish. But those linear terms do not vanish when we 
start from the discrete equation with < 1 (while the quadratic terms do 

not change). In fact, they are divergent in the limit A — > 0. Moreover they are 
"scheme-dependent", in the sense that, even though all cases \i — j\ < 1 lead 
to an equation for (n 2 (x)), this equation is not uniquely obtained; it depends 
on the particular discrete equation that one starts with. We will not write all 
the details here, but we shall give the most divergent terms. We easily find 
that 



d{n 2 (x, x)) 



dY 



div 



4 (n(x)) + O from i = j 

-2 (n(x)) + O (^j from i = j±l. 



(A.15) 



From the previous discussion, we expect that the divergencies will disappear 
when we consider the evolution of normal ordered densities. Indeed, using 
Eq.(A.13) and 



n^(xi,x 2 , x 3 ) = (n(xi)n(x2)n(x 3 )) — [8(xi — x 2 )(n(xi)n(x 3 )} + perm.] 



+ 2 5(a;i - x 2 )S(x 1 - x 3 )(n(x 1 )), 



(A.16) 



and along with Eq.(A.14) we obtain a unique and finite evolution equation of 
the pair density in the elegant form 



<9n (2) (xi,^) 
dY 



= < a 



n^ 2 \xi,x 2 ) + 5(xi — x 2 ) n^\xi) 
- (3 [n ( - 3 \x 1 ,x 1 ,x 2 ) + S(x 1 - x 2 ) n < - 2 \x 1 ,x 1 ) 
+ DnM(x 2 ) ^ i n + {i^ 2 }. 



(A.17) 
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In principle one could follow the same procedure to derive the evolution equa- 
tion for the general K-point function n^(xi, ...,x K ) or understand the form by 
looking at the first two equations, (A. 14) and (A. 17), which we just derived. 
We shall not write the full hierarchy here, but the easiest and most rewarding 
way to derive it is from the corresponding Langevin equation. This is simply 
the generalization of the Langevin equation we saw in Sec. 4 and reads 

= a h(x)-(3h 2 (x)+D + ^2 [a n(x) - (3 n 2 (x)] u(x, Y), (A.18) 

where the noise satisfies (z/(x,F)) = and (u(xi, Y)v(x 2 , Y')) = S(xi — 
x 2 )5(Y — Y'). Notice that h(x) saturates at the value at/ (3. The hierarchy 
generated by Eq.(A.18) is equivalent to the one of the model we studied, pro- 
vided we simply identify (n(xi)...n(x K )) with n^(xi, x K ). One sees that 
in a straightforward way by deriving Eqs.(A.14) and (A. 17) from Eq.(A.18), 
while a formal proof can be found in [38] . 

For our convenience, and in order to make a closer analogy to the QCD 
problem (where T saturates to one), we can make the change of variables 
h — > (a//3)u, Y — ► t/a, D — > aD and v — ► y/a v. Then the Langevin equation 
becomes 



^ = u(x) - u\x) + D ^1 + [u ( x ) - U 2 {x)] ^ t), (A.19) 

with (v(x, t)) = and (v(xi, t)u(x2, t')) = 5{x% — s 2 )5(t — t'). Now u(x) satu- 
rates at the value 1. 
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